# colors for gender
sex_colors =  c(men = "#30D5C8",women = "purple",
                `women-men` = "orange", avg = "blue4")
scale_col_sex = list(
  scale_fill_manual(values = sex_colors), 
  scale_color_manual(values = sex_colors)  
)

sex_colors.b =  c(men = "#30D5C8",women = "purple",
                `Sex difference` = "orange", total = "blue4")
scale_col_sex.b = list(
  scale_fill_manual(values = sex_colors.b), 
  scale_color_manual(values = sex_colors.b)  
)

sex_colors2 =  c("#30D5C8","purple")
scale_col_sex2 = list(
  scale_fill_manual(values = sex_colors2), 
  scale_color_manual(values = sex_colors2)  
)
# colors for drug
drug_colors =  c("blue","grey25","orange")
scale_col_drug = list(
  scale_fill_manual(values = drug_colors), 
  scale_color_manual(values = drug_colors)  
)
# colors for stress condition
stress_colors =  c("#1b9e77","#d95f02")
scale_col_stress = list(
  scale_fill_manual(values = stress_colors), 
  scale_color_manual(values = stress_colors)  
)
gg_no_y_axis = 
  theme(axis.text.y = element_blank(),
        axis.ticks.length.y = unit(0,"mm"),
        panel.grid.major.y = element_blank()) 
gg_no_x_axis = 
  theme(axis.text.x = element_blank(),
        axis.ticks.length.x = unit(0,"mm"),
        panel.grid.major.x = element_blank()) 

negaffect_items = c("distressed","anxious","vulnerable","irritable") 
posaffect_items = c("good", "happy", "confident", "safe", "relaxed")
embselfcn_items = c("embarassed","selfconscious")
sideffect_items = c("dizzy","dry_mouth", "blunted", "nauseous","warm_face","not_myself") 

1 Prepare data

We load the data. The data file was created on NA.

load("data/my_data_Questionnaire.Rdata")
my_data.Q = selector(my_data)
my_data.Q[item_name == "dull", item_name := "blunted"]

load("data/my_data_cortisol.Rdata")
my_data.cortisol = 
  my_data %>% 
  .[, result := result*100*0.0001] %>% 
  .[, log.result := log(result)] %>% 
  .[participant %in% unique(my_data$participant)]

load("data/my_data_SelfAdmin.Rdata")
my_data = selector(my_data)

stage.names = c("BL 1", "BL 2", "BL 3", "state induction", "dose 1" , 
           "reminder 1", "self- admin", "reminder 2", "DM task",
           "cognitive task", "dose 2") %>% gsub(" ","\n",.)
my_data %>% 
  .[, state.condition := ifelse(Stress == "stress", "stress", "control")] %>% 
  .[, Sex := ifelse(gender == "m","men","women")] 
my_data.Q %>% 
  .[, state.condition := ifelse(Stress == "stress", "stress", "control")] %>% 
  .[, Sex := ifelse(gender == "m","men","women")] %>% 
  .[, induced.state := ifelse(Stress == "stress" & stage > 3, "stress", "control")] %>% 
  .[, Stage := factor(stage.names[stage], levels = stage.names)]

Calculate first time at target for self admin:

my_data %>% 
  .[, first.time.target := 1e5 + 0.1] %>% 
  .[, `:=`(first.time.target = 
          get_ftt(effect_obtained,target, time, bias),
          mean.delta_effect_obtained = mean(abs(effect_obtained-target))),
        by = .(participant,session)]

Currently 63 participants (32 women) with data from 4 opioid sessions are included in the analysis.

1.1 Exclusions based on ratings

We exclude people who reported neither an increases in negative emotions nor a decrease in positive emotions after the stress condition.

negaffect_items = c("distressed","anxious","vulnerable","embarassed","irritable") ## add stressed

posaffect_items = c("good", "happy", "confident", "safe", "relaxed")

# maximum change in any of the negative emotion items
my_data.negaffect = 
  my_data.Q %>% 
  .[stage %in% c(3,4) & item_name %in% negaffect_items & state.condition == "stress" & Drug == "oxycodone", 
    .(participant, item_name, response, stage, Sex)] %>% 
  dcast( participant + item_name + Sex ~ stage, value.var = "response") %>% 
  .[, delta := `4` - `3`] %>% 
  .[, .(max.neg_incr = max(delta)), by = .(participant, Sex)]

# maximum change in any of the positive emotion items
my_data.posaffect = 
  my_data.Q %>% 
  .[stage %in% c(3,4) & item_name %in% posaffect_items & state.condition == "stress" & Drug == "oxycodone", 
    .(participant, item_name, response, stage, Sex)] %>% 
  dcast( participant + item_name + Sex ~ stage, value.var = "response") %>% 
  .[, delta := `4` - `3`] %>% 
  .[, .(max.pos_decr = min(delta)), by = .(participant, Sex)]


my_data.posng_change = 
  merge(my_data.negaffect,
        my_data.posaffect,
        by = c("participant","Sex")) %>% 
  setkeyv(c("participant")) %>% 
  .[, stress.filter := ifelse(max.neg_incr <= 0 & max.pos_decr >= 0,1,0)]

my_data.posng_change %>% 
  ggplot(aes(x = max.neg_incr, y = max.pos_decr, color = Sex, label = participant)) + 
  geom_text() + 
  theme(legend.position = "right") +
  scale_col_sex +
  ylab("mac increase negative emotion items") +
  xlab("mac decrease positive emotion items")
Max changes in self-reported negative and positive items from pre to post stress induction.

Figure 1.1: Max changes in self-reported negative and positive items from pre to post stress induction.

0 participants had to be excluded due to insensitivity to the stress induction.

my_data %>% 
  .[my_data.posng_change, stress.filter := stress.filter] 
my_data = my_data[stress.filter == 0]

my_data.Q_keys = key(my_data.Q)
my_data.Q %>% 
  setkeyv(c("participant")) %>% 
  .[my_data.posng_change, stress.filter := stress.filter] %>%  
  setkeyv(my_data.Q_keys)
my_data.Q = my_data.Q[stress.filter == 0]

rm(my_data.posaffect,my_data.posng_change, my_data.negaffect)

Some overview over participants:

  • number participants with self admin data: 63
  • number of missing sessions from self admin data: 0
  • number participants with questionnaire data: 63
  • number of missing sessions from questionnaire data: 1

1.2 Following instructions in the self administration task

We check if participants actually work to obtain the stated desired effect. Participants who end at more than 50 points away from the their declared goal are considered to not have followed the instructions. For these participants the questionnaire data are removed from the analysis.

my_data %>% 
  .[, delta_target := abs(tail(effect_obtained,1) - target[1]), by = .(participant, session)] %>% 
  .[, cheater := delta_target > 50]

tmp = 
  my_data[, .(participant,state.condition,session,Drug, cheater)] %>% 
  unique()

cheaters = 
  tmp %>% 
  .[cheater == TRUE] %>% 
  unique() %>% 
  .[, session := factor(session)] 

cheaters  %>% 
  copy() %>% 
  setnames(c("cheater"),c("Ignored instructions")) %>% 
  .[Drug == "oxycodone"] %>% 
  kable() %>% 
  add_footnote(paste("In",tmp[Drug == "placebo" & cheater == TRUE] %>% nrow(),"of the placebo sessions a participant ignored the instructions.")) %>% 
  kable_styling(full_width = FALSE)
participant state.condition session Drug Ignored instructions
14 control 3 oxycodone TRUE
14 stress 4 oxycodone TRUE
15 control 4 oxycodone TRUE
33 control 3 oxycodone TRUE
33 stress 4 oxycodone TRUE
63 control 3 oxycodone TRUE
63 stress 4 oxycodone TRUE
77 stress 1 oxycodone TRUE
77 control 2 oxycodone TRUE
97 control 1 oxycodone TRUE
120 stress 3 oxycodone TRUE
a In 14 of the placebo sessions a participant ignored the instructions.

In 11 of the () oxycodone sessions the obtained effect was more than 50 points from the target away. Altogether 7 participants contributed such sessions. In 14 of the placebo sessions a participant ignored the instructions.”

tmp %>% 
  .[, .(n = sum(!cheater), p = mean(!cheater)), 
    by = .(Drug, state.condition)] %>%
  .[, N := paste0(n, " (", round(p*100),"%)")] %>% 
  dcast(Drug ~ state.condition, value.var = "N") %>% 
  kable(caption = "Sample size for analyses of drug wanting after removal of participants who did not follow instructions") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic() %>%  
  add_footnote("The number in parentheses indicates the percentage of participant who followed instructions.", notation = "none")
Table 1.1: Sample size for analyses of drug wanting after removal of participants who did not follow instructions
Drug control stress
oxycodone 57 (90%) 58 (92%)
placebo 55 (87%) 57 (90%)
The number in parentheses indicates the percentage of participant who followed instructions.

2 Analysis plan

2.1 Confirmatory analyses

Due to the within-participant manipulation of the independent variables (drug * state) and the non-normal distributions of several of the outcome variables we will use hierarchical Bayesian regression analyses.

In these analyses, random effects will account for repeated measurement of participants (and time whenever measures are repeated within participant). To account for the non-normal distribution of outcome variables (preliminary analysis of pilot data showed some heteroscedasticity of residual from a gaussian / normal model) we will use an item response theory approach to model the data when appropriate.

See the sections Linear regression and Linear regression without sex effects for an analyses of predictions and residuals of linear models that shows that a linear models are not appropriate for rating data or behavioral data from the self administration task. Note that all outcomes are measured on continuations scales, and we use ordinal logistic regressions merely as a flexible modelling approach that allows us to accurately capture patterns in the data. Because the analysed data is generally continuous, it is possible and permissible to calculate expected responses and derived contrasts from the results of an ordinal logistic regression model.

As a preparation for using and ordinal model, we coarsen the data by binning the 0-100 VAS scale into 20 bins of 5 points and the data of the self-administration task, which has a scales of 0-125, into 25 bins or response levels. The motivation for this bin sizes is that it allows to represent questionnaire responses and behavior relatively fine-grained while being coarse enough to allow for relatively easy estimation of threshold parameters of ordinal logistic models.

We use the more flexible adjacent category ordered logistic model because this allows to relax the proportional odds assumptions if needed by estimating category-specific effects. The priors for the thresholds for this model are set as normal distributions with the expected threshold given no experimental effects as mean and2 as the standard deviation. For examples, given four ordered categories with frequencies 10, 20, 40, 30, the expected threshold between categories 1 and 2 is log((10/70)/(20/70)) = -0.3 and the expected threshold between categories 4 and 4 is log((40/70)/(30/70)) = 0.125. Additionally, we use weakly informative priors of N(0,2) for the standard deviation of random effects and N(0,2). We also use weakly informative priors for fix effects of interest. See the section Prior sensitivity analysis in the Main Results for a prior sensitivity analysis for the main outcome.

All models are estimated with the R package brms, which uses the Hamiltonian Markov Chain Monte Carlo Sampler implemented in the programming language Stan to obtain posterior distribution of model parameters. Sampling is generally done with 4 chains, each of which uses 1000 samples for warm up and 1000 post-warmup samples. Warmup samples are discarded and all reported results rely on post-warmup samples. We verify model convergence by checking that Rhat values are below 1.001 and the the effective sample size is above 500 (it is usually above 2000). If 1000 post-warmup samples were not sufficient to achieve 500 effective samples, we re-ran the model with 2000 post-warmup samples.

We report all results by calculating effects (means and 95% credible intervals) on the original scale. We describe effects whose credible intervals do not include 0 as “present” and report in addition to means and credible intervals the posterior probability that an effect is larger than 0 (\(P_{E>0}\)). For ordered logistic regressions we calculate effects as follows: (1) We calculate the average response \(\mu_l\) for each category \(l\) of the coarsened variable from the raw data. (2) Using the posterior distribution from the fitted model, we calculate for each condition \(k\) the level-probabilities \(P_{l,k}\), i.e. the probability with which a response in each category is expected. 1 (3) We calculate the expected response \(R_k\) for each condition as the weighted average of \(mu\) with the level-probabilities \(P_{l,k}\) as weights. (4) We calculate effects of interest as contrasts between the expected responses in different conditions. For instance, stress effects are generally calculated as changes in stress response from pre to post state induction.

This regression analysis also allows adjusting for covariates such as session, drug and state induction order, gender, age, BMI, body weight, time of day, experimenter/medical staff identity, drug side effects.

The our analyses adjust for session number, whereby we model the categorical sessions variable hierarchically, in order to regularize session effects. Given that preliminary analyses showed strong sex effects, we also modeled main and interaction effects for sex. For these analysis we generally calculated sex-specific and over-all (averaged over sex) effects as well as sex differences. However, consistent with the pre-registration, we also estimate a model without an interaction with sex for the main outcome from the self administration task. Where relevant, we model repeated measurement of the same item with random effects.

An initial check of the data also revealed substantial effect heterogeneity, i.e. there was great variation in the response to stress and drug manipulations. We took this into account by explicitly modeling it as random treatment effects. Indeed, model comparison of models with and without random treatment effects showed that the latter has a better out of sample predictive value.

Main predictors: Drug condition, stress condition (Note that some of the hypotheses listed below are tested within one of the drug or stress conditions only).

Please see the supplementary materials for complete specifications of all analyses and code that can be used to re-produce the results.

3 Manipulation check

3.1 Stress

3.1.1 Linear regression

my_data.Q.stressed = select_data("stressed")
my_stages = c(1:5)

my_data.Q.stressed %>% plot_stages()
Stress over time. Lines are mean responses. Confidence intervals are within-confidence intervals.

Figure 3.1: Stress over time. Lines are mean responses. Confidence intervals are within-confidence intervals.

We measure stress as the change in stress rating from pre to post induction. The following plot shows this:

the_data = make_the_data(my_data.Q.stressed)[, outcome := "stress"]
attr(the_data,"N")
Table 3.1: Number of particicipants per condition by condition and stage
Drug state.condition pre.induction post.induction
oxycodone control 62 62
oxycodone stress 63 63
placebo control 63 63
placebo stress 63 63
the_data %>% plot_pre_post() + scale_col_stress
Change in stress ratings from pre to post state induction. Each line represents an individual. Stress is reflected in the slope of the lines.

Figure 3.2: Change in stress ratings from pre to post state induction. Each line represents an individual. Stress is reflected in the slope of the lines.

We can further summarize this information by showing the distribution of the slopes in the previous figure:

p = plot_stresseff_raw(the_data)
raw_contrasts.stress = 
  p$data %>% 
  .[, .(participant,Sex,Drug, stress_effect)] %>% 
  .[, Outcome := "Stressed (VAS)"]
p
Distribution of change in stress ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Figure 3.3: Distribution of change in stress ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

The next plot shows the individual data points, where each data point is calculated as a pre-post difference, together with mean and standard errors.

plot_stresseff_raw2(the_data)
Mean and CIs of change in stress ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Figure 3.4: Mean and CIs of change in stress ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

To estimate the effect of stress induction on stress ratings, we use a hierarchical ordered logistic model regression. The dependent variables are pre and post stress induction ratings and the independent variable the state induction. We model the responses as main effects of Stress and Sex plus the state-induction X sex interaction, where the interaction term measures if the effect of stress induction differed between males and females:

stress_formula = response.c ~ induced.state + Sex + stage:induced.state:Sex + (1|session) + (1|participant)

stress_formula = response.c ~ induced.state + Sex + stage:induced.state:Sex + (1|session) + (1|participant) 

We estimate a standard hierarchical linear regression with Gaussian errors.

ffit = lme4::lmer(response ~ induced.state + Sex + stage:induced.state:Sex + (1|session) + (1|participant) , data = the_data)

To check the model, we first investigate model predictions. The following figure shows that the model predictions, and even the best linear predictor predict impossible negative values.

error_sd = 
  VarCorr(ffit) %>% attr("sc")
p1 = 
  data.table(yhat = predict(ffit)) %>% 
  .[, lower := yhat - error_sd] %>% 
  .[, upper := yhat + error_sd] %>% 
  .[, state.condition := the_data$state.condition] %>% 
  .[order(yhat)] %>% 
  .[, j := 1:.N] %>% 
  ggplot(aes(x = j, y = yhat, ymin = lower, ymax = upper,
             color = state.condition, fill = state.condition)) + 
   geom_hline(yintercept = 0, color = "red") +
  geom_linerange(alpha = .5) +
  geom_point(size = .4) +
  xlab("") +
  ylab("fitted + error") + 
  scale_col_stress
p2 = 
  rbind(the_data[, .(response)][, data := "observed"],
           data.table(response = predict(ffit), data = "predicted")) %>% 
  ggplot(aes(x = response, fill = data)) + 
  geom_histogram(position = "identity", alpha = .5)
  
(p1 | p2) + plot_layout(widths = c(2,1))
Predictions of the linear model +/- 1sd of the normal error distribution.

Figure 3.5: Predictions of the linear model +/- 1sd of the normal error distribution.

Consistently, the following figure shows that the residuals are correlated with the outcome variable and not normally distributed.

res = data.table(residuals = residuals(ffit), fitted = predict(ffit))
p1 = res %>% 
  ggplot(aes(x = fitted, y = residuals)) + 
  geom_point() + 
  geom_smooth()
p2 = 
  res %>% 
  ggplot(aes(x = residuals)) + 
  geom_histogram(aes(y =..density..),
                 breaks = seq(-60, 60, by = 5), 
                 colour = "black", 
                 fill = "white") +
  stat_function(fun = dnorm,
                args = list(mean = mean(res$residuals),
                            sd = sd(res$residuals)))

p1 | p2
Residuals from a linear regression model. Left: Predicted values and residuals are not independent. Right: Residuals are not normally distributed.

Figure 3.6: Residuals from a linear regression model. Left: Predicted values and residuals are not independent. Right: Residuals are not normally distributed.

rm(error_sd,ffit,p1,p2,res)

In sum, a linear regression model with Gaussian errors is not a good model for the rating scale data. We therefore use ordinal regression models for coarsened data. Specifically, we coarsen the questionnaire data into 20 equally sized bins and perform Bayesian ordinal logistic regressions on the coarsened data. To report results, we back-calculate effects (including credible intervals) of independent variables on the original 0-100 rating scale.

3.1.2 Ordinal logistic regression

We coarsen the outcome variable in preparation of the ordinal regression model

the_data = the_data %>% 
  coarsen(var = "response", bins = 20, range = c(0,100))

Here we run the Bayesian ordinal logistic regression for the first stress induction:

fn = "brmsP1/stressed_acat.Rdata"
if (file.exists(fn)) {
  load(fn)
} else {
  fit = fit_ordinal_model(the_data, stress_formula, fn, family = acat())
}
Click to show sample diagnostics
Table 3.2: Summary of model paramters
parameter mean sd q5 q95 rhat ess_bulk ess_tail
Threshold[1] 1.3 1.13 -0.6 3.1 1.001 1729 2599
Threshold[2] 0.5 1.13 -1.4 2.3 1.001 1673 2187
Threshold[3] 0.4 1.15 -1.5 2.3 1.001 1765 2496
Threshold[4] 1.1 1.17 -0.9 3.0 1.002 1674 2509
Threshold[5] 0.3 1.18 -1.7 2.2 1.001 1823 2714
Threshold[6] 0.8 1.20 -1.2 2.7 1.001 1817 2730
Threshold[7] 0.6 1.23 -1.5 2.6 1.000 1809 2527
Threshold[8] 0.7 1.23 -1.3 2.7 1.001 1880 2510
Threshold[9] 0.4 1.22 -1.6 2.4 1.001 1831 2586
Threshold[10] 0.5 1.22 -1.4 2.5 1.001 1832 2726
Threshold[11] 0.8 1.22 -1.2 2.8 1.001 1797 2584
Threshold[12] 1.1 1.22 -1.0 3.1 1.001 1906 2855
Threshold[13] 0.8 1.25 -1.3 2.9 1.002 1903 2662
Threshold[14] 0.5 1.24 -1.6 2.6 1.001 1890 2659
Threshold[15] 0.8 1.22 -1.3 2.8 1.002 1826 2567
Threshold[16] 1.5 1.26 -0.5 3.6 1.001 1957 2567
Threshold[17] 1.2 1.29 -0.9 3.3 1.001 1998 2616
Threshold[18] 1.0 1.31 -1.2 3.2 1.001 1944 2541
Threshold[19] 0.7 1.30 -1.5 2.8 1.001 1975 2624
b_induced.statestress 0.4 1.33 -1.8 2.6 1.001 2033 2460
b_Sexwoman 0.1 1.23 -2.0 2.0 1.002 1956 2470
b_induced.statecontrol:Sexman:stagepre.induction -0.2 1.11 -2.0 1.6 1.001 1695 2674
b_induced.statestress:Sexman:stagepre.induction 0.0 2.01 -3.4 3.3 1.000 6573 2797
b_induced.statecontrol:Sexwoman:stagepre.induction 0.0 1.14 -1.9 1.8 1.003 1612 2288
b_induced.statestress:Sexwoman:stagepre.induction 0.0 2.06 -3.4 3.4 1.000 5690 2653
b_induced.statecontrol:Sexman:stagepost.induction -0.2 1.11 -2.0 1.6 1.001 1702 2693
b_induced.statestress:Sexman:stagepost.induction 0.1 1.34 -2.1 2.4 1.001 2249 2742
b_induced.statecontrol:Sexwoman:stagepost.induction -0.1 1.14 -1.9 1.8 1.002 1611 2316
b_induced.statestress:Sexwoman:stagepost.induction 0.3 1.37 -1.9 2.5 1.003 2147 2807
sd_participant__Intercept 0.3 0.05 0.2 0.4 1.002 1429 2363
sd_session__Intercept 0.0 0.07 0.0 0.2 1.005 1290 1338
model: response.c ~ induced.state + Sex + stage:induced.state:Sex + (1 | session) + (1 | participant)
b_ = fix effect parameters
sd_ standard deviation of hierarchical parameters
ess = effective sample size


The next figure shows that the model describes the observed data well.

#pp_check(fit, ndraws = 4000, type = "bars")

Due to non-collapsibility of OR we cannot simply interpret the interaction effect. Instead, We calculate the contrasts manually from posterior predictions:

tmp = plot_prepost_contrast(fit, the_data,"stress", my_contrast = "stress")
tmp[[1]]
Posterior distributions of estimated stress induction effects on self reported stress on the 0-100 VAS scale in women, men, and sex difference (women-men). Density plots show posterior distributions. Dots and lines identify means, and 50 & 95% credible intervals.

Figure 3.7: Posterior distributions of estimated stress induction effects on self reported stress on the 0-100 VAS scale in women, men, and sex difference (women-men). Density plots show posterior distributions. Dots and lines identify means, and 50 & 95% credible intervals.

tmp$table
Table 3.3: Modeled post minus pre state induction on stress ratings by condition (control, stress) and the stress effect by sex.
Sex control stress stress_effect
men 0 [-2, 3]; 29 [24, 34]; 29 [24, 34];
women -0 [-3, 3]; 47 [42, 52]; 48 [42, 53];
total -0 [-2, 2]; 38 [34, 42]; 38 [34, 42];
Sex difference -1 [-4, 3]; 18 [10, 25]; 18 [11, 26];
post = tmp[[2]]
summary_table.stress = tmp$summary_table[, Outcome := "Stressed (VAS)"]

The estimated effect of stress from the state induction is NaN (CI=[NA, NA]; BFNA; P(S>C)NA). The effect is 48 (CI=[42, 53]; BF>10^5; P(S>C)>0.999) in females, 29 (CI=[24, 34]; BF>10^5; P(S>C)>0.999) in males and the sex difference is NaN (CI=[NA, NA]; BFNA; P(S>C)NA).

3.2 Negative mood

We repeat the same analysis as for stress and also take into account that there are multiple items measuring negative mood. These items are: distressed, anxious, vulnerable, embarassed, irritable

my_data.Q.neg_mood = select_data(negaffect_items)
my_data.Q.neg_mood %>% plot_stages()
Negative mood over time. Lines are mean responses.

Figure 3.8: Negative mood over time. Lines are mean responses.

the_data = make_the_data(my_data.Q.neg_mood)[, outcome := "neg_mood"]
attr(the_data,"N")
Table 3.4: Number of particicipants per condition by condition and stage
Drug state.condition pre.induction post.induction
oxycodone control 62 62
oxycodone stress 63 63
placebo control 63 63
placebo stress 63 63
the_data %>% plot_pre_post() + scale_col_stress
Change in negative mood ratings from pre to post state induction. Each line represents an individual. Stress is reflected in the slope of the lines.

Figure 3.9: Change in negative mood ratings from pre to post state induction. Each line represents an individual. Stress is reflected in the slope of the lines.

We can further summarize this information by showing the distribution of the slopes in the previous figure:

plot_stresseff_raw(the_data)
Distribution on change in negative mood ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Figure 3.10: Distribution on change in negative mood ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

plot_stresseff_raw2(the_data)
Mean and CIs of change in negative mood ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Figure 3.11: Mean and CIs of change in negative mood ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

We extend the analysis model to account for multiple items.

stress_formula = response.c ~ stage:induced.state:Sex + session + (1|participant) + (1 | item_name)

stress_formula = response.c ~ stage:induced.state:Sex + (1|session) + (1|participant) + (1 | item_name)

Here we run the Bayesian ordinal logistic regression:

fn = "brmsP1/neg_mood_acat.Rdata"
the_data = the_data %>% 
  coarsen(var = "response", bins = 20, range = c(0,100))
if (file.exists(fn)) {
  load(fn)
} else {
  fit = fit_ordinal_model(the_data, stress_formula, fn, family = acat())
}
#pp_check(fit, ndraws = 2000, type = "bars")
Click to show sample diagnostics
Table 3.5: Summary of model paramters
parameter mean sd q5 q95 rhat ess_bulk ess_tail
Threshold[1] 1.2 0.84 -0.2 2.5 1.003 1013 1426
Threshold[2] 0.8 0.84 -0.6 2.1 1.003 1032 1445
Threshold[3] 0.5 0.85 -0.9 1.9 1.003 1052 1352
Threshold[4] 0.6 0.86 -0.9 1.9 1.003 1045 1482
Threshold[5] 0.4 0.87 -1.0 1.8 1.003 1100 1510
Threshold[6] 0.2 0.88 -1.1 1.6 1.003 1084 1523
Threshold[7] 0.9 0.89 -0.6 2.3 1.003 1099 1567
Threshold[8] 0.5 0.92 -1.0 2.0 1.004 1134 1653
Threshold[9] -0.3 0.90 -1.8 1.1 1.003 1150 1534
Threshold[10] 0.8 0.89 -0.7 2.2 1.003 1144 1593
Threshold[11] 1.0 0.92 -0.4 2.5 1.002 1151 1868
Threshold[12] 0.0 0.92 -1.5 1.5 1.003 1163 1557
Threshold[13] 0.4 0.90 -1.1 1.9 1.003 1150 1621
Threshold[14] 1.3 0.92 -0.2 2.8 1.003 1148 1795
Threshold[15] 0.7 0.95 -0.9 2.3 1.003 1215 1731
Threshold[16] 0.5 0.96 -1.1 2.1 1.004 1265 1964
Threshold[17] 2.4 1.11 0.6 4.2 1.001 1541 2185
Threshold[18] 1.5 1.33 -0.7 3.7 1.001 1815 2153
Threshold[19] -1.8 1.25 -4.0 0.1 1.002 1672 2282
b_stagepre.induction:induced.statecontrol:Sexman -0.3 0.82 -1.7 1.0 1.003 1060 1380
b_stagepost.induction:induced.statecontrol:Sexman -0.4 0.82 -1.8 0.9 1.002 1075 1408
b_stagepre.induction:induced.statestress:Sexman 0.0 1.98 -3.3 3.2 1.001 6671 2940
b_stagepost.induction:induced.statestress:Sexman 0.2 0.82 -1.1 1.5 1.003 1060 1341
b_stagepre.induction:induced.statecontrol:Sexwoman -0.1 0.82 -1.4 1.2 1.003 1062 1376
b_stagepost.induction:induced.statecontrol:Sexwoman -0.2 0.82 -1.5 1.1 1.003 1060 1397
b_stagepre.induction:induced.statestress:Sexwoman 0.0 2.03 -3.3 3.4 1.001 6800 2443
b_stagepost.induction:induced.statestress:Sexwoman 0.4 0.82 -0.9 1.7 1.003 1062 1454
sd_item_name__Intercept 0.3 0.22 0.1 0.7 1.002 1180 1233
sd_participant__Intercept 0.3 0.04 0.2 0.3 1.002 1208 1921
sd_session__Intercept 0.1 0.11 0.0 0.3 1.005 1064 1478
model: response.c ~ stage:induced.state:Sex + (1 | session) + (1 | participant) + (1 | item_name)
b_ = fix effect parameters
sd_ standard deviation of hierarchical parameters
ess = effective sample size


We calculate the contrasts manually from posterior predictions:

tmp = plot_prepost_contrast(fit, the_data,"neg_mood", my_contrast = "stress")
tmp[[1]]
Posterior distributions of estimated stress induction effects on self reported negative mood on the 0-100 VAS scale in women, men, and sex difference (women-men). Density plots show posterior distributions. Dots and lines identify means, and 50 & 95% credible intervals.

Figure 3.12: Posterior distributions of estimated stress induction effects on self reported negative mood on the 0-100 VAS scale in women, men, and sex difference (women-men). Density plots show posterior distributions. Dots and lines identify means, and 50 & 95% credible intervals.

tmp$table
Table 3.6: Modeled post minus pre state induction on neg_mood ratings by condition (control, stress) and the stress effect by sex.
Sex control stress stress_effect
men -0 [-1, 0]; 13 [11, 15]; 13 [11, 16];
women -1 [-2, 1]; 25 [22, 28]; 26 [23, 28];
total -1 [-1, 0]; 19 [17, 21]; 20 [18, 21];
Sex difference -0 [-2, 1]; 12 [9, 15]; 12 [9, 16];
post = tmp[[2]]
summary_table.stress = 
  rbind(summary_table.stress,
        tmp$summary_table[, Outcome := "Negative Mood"])

The estimated effect of stress from the state induction is NaN (CI=[NA, NA]; BFNA; P(S>C)NA). The effect is 26 (CI=[23, 28]; BF>10^5; P(S>C)>0.999) in females, 13 (CI=[11, 16]; BF>10^5; P(S>C)>0.999) in males and the sex difference is NaN (CI=[NA, NA]; BFNA; P(S>C)NA).

3.3 Positive mood

We repeat the same analysis as for stress and also take into account that there are multilple items measuring negative mood. These items are: distressed, anxious, vulnerable, embarassed, irritable

my_data.Q.pos_mood = select_data(posaffect_items)
my_data.Q.pos_mood %>% plot_stages()
Positive mood over time. Lines are mean responses.

Figure 3.13: Positive mood over time. Lines are mean responses.

the_data = make_the_data(my_data.Q.pos_mood)[, outcome := "pos_mood"]
attr(the_data,"N")
Table 3.7: Number of particicipants per condition by condition and stage
Drug state.condition pre.induction post.induction
oxycodone control 62 62
oxycodone stress 63 63
placebo control 63 63
placebo stress 63 63
the_data %>% plot_pre_post() + scale_col_stress
Change in positive mood ratings from pre to post state induction. Each line represents an individual. The stress induction effect is reflected in the slope of the lines.

Figure 3.14: Change in positive mood ratings from pre to post state induction. Each line represents an individual. The stress induction effect is reflected in the slope of the lines.

We can further summarize this information by showing the distribution of the slopes in the previous figure:

plot_stresseff_raw(the_data)
Distribution on change in positive mood ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Figure 3.15: Distribution on change in positive mood ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

plot_stresseff_raw2(the_data)
Mean and CIs of change in positive mood ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Figure 3.16: Mean and CIs of change in positive mood ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Here we run the Bayesian ordinal logistic regression:

fn = "brmsP1/pos_mood_acat.Rdata"
the_data = the_data %>% 
  coarsen(var = "response", bins = 20, range = c(0,100))
if (file.exists(fn)) {
  load(fn)
} else {
  fit = fit_ordinal_model(the_data, stress_formula, fn, family = acat())
}
#pp_check(fit, ndraws = 2000, type = "bars")
Click to show sample diagnostics
Table 3.8: Summary of model paramters
parameter mean sd q5 q95 rhat ess_bulk ess_tail
Threshold[1] 0.6 0.89 -0.8 2.1 1.010 659 1230
Threshold[2] -0.5 0.87 -1.9 1.0 1.007 687 1079
Threshold[3] -0.9 0.89 -2.3 0.6 1.009 681 1263
Threshold[4] 0.0 0.87 -1.4 1.5 1.009 662 1224
Threshold[5] -0.5 0.87 -1.9 1.0 1.008 676 945
Threshold[6] -0.5 0.88 -1.8 1.0 1.009 672 1219
Threshold[7] -0.3 0.86 -1.6 1.2 1.010 645 950
Threshold[8] -0.1 0.87 -1.5 1.3 1.008 670 1101
Threshold[9] -0.9 0.86 -2.2 0.6 1.009 645 1094
Threshold[10] 0.0 0.86 -1.3 1.4 1.009 651 1201
Threshold[11] 0.0 0.86 -1.4 1.4 1.009 651 1155
Threshold[12] -0.3 0.85 -1.6 1.2 1.009 645 1059
Threshold[13] 0.1 0.86 -1.3 1.5 1.009 643 1071
Threshold[14] -0.1 0.86 -1.4 1.4 1.009 649 1105
Threshold[15] 0.1 0.85 -1.2 1.5 1.009 643 968
Threshold[16] 0.2 0.85 -1.1 1.7 1.010 637 1203
Threshold[17] 0.2 0.85 -1.2 1.6 1.009 647 1125
Threshold[18] 0.7 0.86 -0.6 2.1 1.009 654 1067
Threshold[19] -0.7 0.86 -2.1 0.7 1.009 640 1150
b_stagepre.induction:induced.statecontrol:Sexman 0.2 0.84 -1.1 1.6 1.011 645 1059
b_stagepost.induction:induced.statecontrol:Sexman 0.2 0.84 -1.1 1.6 1.011 645 1021
b_stagepre.induction:induced.statestress:Sexman 0.0 2.03 -3.4 3.3 1.002 3811 2923
b_stagepost.induction:induced.statestress:Sexman -0.1 0.84 -1.4 1.3 1.011 644 1088
b_stagepre.induction:induced.statecontrol:Sexwoman 0.1 0.84 -1.2 1.5 1.010 650 1136
b_stagepost.induction:induced.statecontrol:Sexwoman 0.1 0.84 -1.2 1.5 1.010 649 1105
b_stagepre.induction:induced.statestress:Sexwoman -0.1 1.96 -3.3 3.2 1.000 3766 2717
b_stagepost.induction:induced.statestress:Sexwoman -0.3 0.84 -1.6 1.1 1.010 650 1144
sd_item_name__Intercept 0.2 0.10 0.1 0.3 1.004 820 1189
sd_participant__Intercept 0.2 0.02 0.2 0.2 1.008 548 1200
sd_session__Intercept 0.0 0.06 0.0 0.1 1.001 768 1024
model: response.c ~ stage:induced.state:Sex + (1 | session) + (1 | participant) + (1 | item_name)
b_ = fix effect parameters
sd_ standard deviation of hierarchical parameters
ess = effective sample size


We calculate the contrasts manually from posterior predictions:

tmp = plot_prepost_contrast(fit, the_data, "pos_mood", my_contrast = "stress")
tmp[[1]]
Posterior distributions of estimated stress induction effects on self reported positive mood on the 0-100 VAS scale in women, men, and sex difference (women-men). Density plots show posterior distributions. Dots and lines identify means, and 50 & 95% credible intervals.

Figure 3.17: Posterior distributions of estimated stress induction effects on self reported positive mood on the 0-100 VAS scale in women, men, and sex difference (women-men). Density plots show posterior distributions. Dots and lines identify means, and 50 & 95% credible intervals.

tmp$table
Table 3.9: Modeled post minus pre state induction on pos_mood ratings by condition (control, stress) and the stress effect by sex.
Sex control stress stress_effect
men 1 [-1, 3]; -19 [-21, -16]; -20 [-23, -17];
women 2 [-1, 4]; -33 [-36, -30]; -35 [-38, -32];
total 1 [-0, 3]; -26 [-28, -24]; -27 [-29, -25];
Sex difference 1 [-3, 4]; -14 [-18, -11]; -15 [-19, -11];
post = tmp[[2]]
summary_table.stress = 
  rbind(summary_table.stress,
        tmp$summary_table[, Outcome := "Positive Mood"])

The estimated effect of stress from the state induction is NaN (CI=[NA, NA]; BFNA; P(S>C)NA). The effect is -35 (CI=[-38, -32]; BF<10^-5; P(S>C)<0.001) in females, -20 (CI=[-23, -17]; BF<10^-5; P(S>C)<0.001) in males and the sex difference is NaN (CI=[NA, NA]; BFNA; P(S>C)NA).

3.4 Heart rate

See the script “HRV.R” for code for pre-processing of the raw heart rate data.

In the data.table HR the column time indicates seconds from the middle of the stress induction and BPS indicates heart beats per second. We down sample to measure heart rate only every 15s.

load("data/HRV/HRV.Rdata")
HR = 
  HR %>% 
  .[, participant := as.numeric(gsub("E","",Participant))] %>% 
  .[participant %in% unique(my_data$participant)] %>% 
  .[, state.condition := ifelse(Stress == "Stress", "stress", "control")] %>% 
  .[, Sex := ifelse(Gender == "Male","men","women")] %>% 
  .[, minute := (time/60)]  

tHR = 
  HR %>% 
  .[, time.c := cut(time,
                    breaks = seq(-1200,1200,by = 20),
                    include.lowest = TRUE)] %>% 
  .[, .(time = mean(time), BPS = mean(BPS)),
    by = .(Participant, Session,Sex, Drug, Condition, state.condition, time.c)] %>% 
  .[, time.c := NULL]
tmp = 
  HR %>% 
  .[abs(minute) < 3.5] %>% # stress induction period
  .[, .(mBPS = mean(BPS)), by = .(participant, Sex, Drug, state.condition)] %>% 
  dcast(participant + Sex + Drug ~ state.condition, value.var = "mBPS") %>% 
  .[, stress_effect := stress-control] %>% 
  .[, .(participant, Sex, Drug, stress_effect)] %>% 
  .[, Outcome := "Heart rate, BPM"]

raw_contrasts.stress = 
  rbind(
    raw_contrasts.stress,
    tmp
  )
rm(tmp)

To estimate the effect of stress we fit a Bayesian spline model and compare the estimated heart rate after removing session effects.

Before the model can be estimated, we prepare the date by

  • re-scaling time and heart rate (BPS) to facilitate parameter estimation.
  • creating a group variable grp by crossing the Stress, Drug, and Sex variable. This will allow us to fit splines for these “groups” and makes it possible to estimate stress effects stratified by drug condition and sex.
  • creating a auto-regression group variable AR.grp by crossing the variables participant and session. This variable will be used to group observations for the estimation of autoregressive error terms.
m.BPS = mean(tHR$BPS)
sd.BPS = mean(tHR$BPS)
tHR %>% 
  .[, sBPS := (BPS-m.BPS)/sd.BPS] %>% 
  .[, time := time/1000] 

tHR %>% 
  .[, grp := paste(c(as.character(state.condition),Drug,Sex), collapse = ":"), by = 1:nrow(tHR)] %>% 
  .[, AR.grp := paste(c(Participant,Session), collapse = ":"), by = 1:nrow(tHR)]

The brms model is then specified as:

HRV.formula = sBPS ~ 
  s(time, by = grp) + # smooth effects stratified by sex, drug, and stress condition
  (1|Session) +
  (1|Participant) + # participant and session level random effects
  ar(time = time, gr = AR.grp, p = 1, cov = FALSE) # autoregressive error terms

This model includes an adjustment for session, nested in participants.

Next we estimate the model:

fn = "brmsP1/HRV_induction_2.Rdata"

if (file.exists(fn)) {
  load(fn)
} else {
  bfit =  
    brm(HRV.formula,
        family = gaussian(),
        backend = "cmdstanr",
        iter = 1000, 
        data = tHR,
        prior = prior(normal(0,2), class = "b"))
  save(bfit, file = fn)
}
Click to show sample diagnostics
Table 3.10: Summary of model paramters
parameter mean sd q5 q95 rhat ess_bulk ess_tail
b_Intercept 0.0 0.02 0.0 0.0 1.027 168 366
bs_stime:grpcontrol:Oxycodone:man_1 0.0 0.40 -0.6 0.7 1.004 1211 1347
bs_stime:grpcontrol:Oxycodone:woman_1 -0.3 0.38 -0.9 0.3 1.003 1092 1259
bs_stime:grpcontrol:Placebo:man_1 0.0 0.41 -0.6 0.7 1.001 859 1209
bs_stime:grpcontrol:Placebo:woman_1 0.4 0.35 -0.2 0.9 1.002 1321 1267
bs_stime:grpstress:Oxycodone:man_1 -0.1 0.42 -0.8 0.6 1.004 1259 1202
bs_stime:grpstress:Oxycodone:woman_1 0.2 0.41 -0.5 0.8 1.000 1341 1391
bs_stime:grpstress:Placebo:man_1 0.8 0.38 0.2 1.5 1.001 1060 1395
bs_stime:grpstress:Placebo:woman_1 0.0 0.40 -0.6 0.7 1.001 1404 1268
sd_Participant__Intercept 0.1 0.01 0.1 0.1 1.029 227 411
sd_Session__Intercept 0.0 0.02 0.0 0.0 1.003 447 488
sigma 0.1 0.00 0.1 0.1 1.007 5462 1202
model: sBPS ~ s(time, by = grp) + (1 | Session) + (1 | Participant) + ar(time = time, gr = AR.grp, p = 1, cov = FALSE)
b_ = fix effect parameters
sd_ standard deviation of hierarchical parameters
ess = effective sample size


Here are the expected heart rates split by state and drug conditions and sex:

fn = "brmsP1/HRV_induction_2_pp.Rdata"

# new data for model predictions, where all participants have the same times
new.data = unique(tHR[, .(grp,Participant, Session, AR.grp)])
times = seq(min(bfit$data$time), max(bfit$data$time), length.out = 51)
new.data = 
  do.call(rbind,
          lapply(1:nrow(new.data), 
                 function(x) cbind(new.data[x],time = times)))

# generate posterior predictions
if (file.exists(fn)) {
  load(fn)
} else {
  pp = posterior_epred(bfit, newdata = new.data, cores = 4)
  
  # merge predictions with  drug, stress, sex and time info
  pdata = cbind(new.data, t(pp)) %>% 
    melt(id.vars = names(new.data), variable.name = "iter") %>% 
    .[, iter := gsub("V","",iter)]
  
  # average over participants with the same drug, stress, sex values
  pdata = pdata %>% 
    .[, .(value = mean(value)), by = .(grp,iter, time)]
  
  # prepare data for plotting
  pdata %>% 
    .[, c("Stress","Drug","Gender") := tstrsplit(grp,":")] %>% 
    .[, state.condition := factor(Stress, levels = c("stress","control"))] %>% 
    .[, Drug := factor(Drug)] %>% 
    .[, Sex := factor(Gender)] %>% 
    .[, iter := factor(iter)] %>% 
    # bring outcome variable back to original scale
    .[, value := value*sd.BPS+m.BPS] 
  
  save(new.data, pp , pdata, file = fn)
}


# create a baseline corrected version of the outcome,
# where the drug/stress/sex "groups" have the same
# average heart rate in the minutes -20 to -10
pdata = 
  pdata %>% 
  .[, minute := (time*1000/60)] %>% 
  .[, bl := mean(value[minute < -10]), by = .(state.condition, Drug, Gender, iter)] %>% 
  .[, value.blc := value-bl+m.BPS] %>% 
  .[, sc.i := paste(state.condition,iter)]

pdata[, Sex := gsub("man","men",Sex)]

pdata = rbind(
  pdata,
  avg = pdata %>% 
    .[,
      .(value = fmean(value), bl = fmean(bl), value.blc = fmean(value.blc)),
      by = .(Stress,Drug,state.condition,time,minute,iter,sc.i)] %>% 
    .[, Sex := "avg"],
  fill = TRUE)
# plot the expected heart rates from the first 100 iterations
pdata %>% 
  .[iter %in% 1:100] %>% 
  .[, alpha := 1] %>% 
  ggplot(aes(x = minute, y = value.blc, color = state.condition, group = sc.i, alpha = alpha)) + 
  geom_vline(xintercept = c(-3.5,3.5)) +
  geom_line() + 
  scale_alpha(range = c(0,.1),guide="none") + 
  facet_grid(Sex~Drug) +
  xlab("Time in experiment") + 
  ylab("Beats per minute") + 
  scale_x_continuous(breaks = c(-20,0,20), labels = c("BL3","state induction","dose1")) + 
  scale_col_stress + 
  theme(legend.position = "top")
Estimated heart rate as a function of drug and state condition and sex. The black vertical lines indicate the start and the end of the stress induction period.

Figure 3.18: Estimated heart rate as a function of drug and state condition and sex. The black vertical lines indicate the start and the end of the stress induction period.

Next we calculate the stress effects by comparing heart rates during the induction period:

test.data = 
  pdata %>% 
  .[abs(minute)<3.5, # limit to induction period
    .(value.blc = mean(value.blc)), # calculate average of baseline corrected HR
    by = .(state.condition, Drug, Sex, iter)] %>% 
  # transform to long format
  melt(id.vars = c("state.condition", "Drug", "Sex", "iter"), variable.name = "outcome") %>%   # transform to wide format, so that one column is stress and the other control
  dcast(iter + Drug + Sex ~ state.condition, value.var = "value") %>% 
  # calculate stress effect
  .[, Stress_effect := stress - control]

sexdiff = 
  test.data %>% 
  .[Sex != "avg"] %>% 
  dcast(iter + Drug ~ Sex, value.var = "Stress_effect") %>% 
  .[, Stress_effect := women-men] %>% 
  .[, Sex := "women-men"]

test.data = rbind(
  test.data, sexdiff,
  fill = TRUE
)

test.data %>% 
  ggplot(aes(x = Stress_effect, fill = Sex)) + 
  stat_halfeye(alpha = .5) + 
  facet_grid(Drug~Sex) + 
  geom_vline(xintercept = 0) +
  scale_col_sex + 
  ylab("") + gg_no_y_axis
Posterior distributions of the effect of stress induction on heart rate stratified by drug condition and sex. Density plots show posterior distributions. Dots and lines identify means, and 50 & 95% credible intervals.

Figure 3.19: Posterior distributions of the effect of stress induction on heart rate stratified by drug condition and sex. Density plots show posterior distributions. Dots and lines identify means, and 50 & 95% credible intervals.

stress.effect = test.data[, .(Stress_effect = mean(Stress_effect)), by = .(Sex,iter)]
# update table with stress effects
old.nms = c("men", "women", "avg", "women-men")
new.nms = c("men", "women", "total", "Sex difference")
tmp = stress.effect[, .(get_mci(Stress_effect, get.P = ifelse(grepl("-",Sex),TRUE,FALSE))), by = .(Sex)] %>% 
  dcast(.~Sex) %>% 
  setnames(old.nms,new.nms) %>% 
  .[,new.nms,with = FALSE]
summary_table.stress = rbind(
  summary_table.stress,
  cbind(Outcome = "Heart rate, BPM",tmp)
)
rm(tmp)

On average, stress increases heart rate during the state induction by 12.4 (CI=[9.9, 14.9]; BF>10^5; P(S>C)>0.999) by beats per minutes. The effect in females was 15.1 (CI=[11.6, 18.5]; BF>10^5; P(S>C)>0.999) and 9.8 (CI=[6.4, 13.4]; BF>10^5; P(S>C)>0.999) in males the effect of sex on the stress effect is 5.2 (CI=[0.3, 10.5]; BF=47.8; P(S>C)=0.98).

Lastly, we formally test if the stress effects differs between drug conditions or between sexes.

contrast.data = 
  rbind(
  test.data %>% 
  dcast(iter + Drug ~ Sex, value.var = "Stress_effect") %>% 
  .[, estimate := women-men] %>% 
  .[, .(estimate = mean(estimate)), by = .(iter)] %>% 
  .[, contrast := "Stress:Sex (women-men)"],
  test.data %>% 
  dcast(iter + Sex ~ Drug, value.var = "Stress_effect") %>% 
  .[, estimate := Oxycodone-Placebo] %>% 
  .[, .(estimate = mean(estimate)), by = .(iter)] %>% 
  .[, contrast := "Stress:Drug (Oxycodone-Placebo)"])

contrast.data %>% 
  ggplot(aes(x = estimate)) + 
  geom_vline(xintercept = 0) +
  stat_halfeye(alpha = .5) + 
  facet_wrap(~contrast, scales = "free_x") +
  xlab("Estimated interaction effect") + 
  ylab("") + gg_no_y_axis
Interaction effects. Left: No clear differences in the effect of state induction between drug conditions. Right: State induction increases heart rate in females more than in males.

Figure 3.20: Interaction effects. Left: No clear differences in the effect of state induction between drug conditions. Right: State induction increases heart rate in females more than in males.

The stress induced heart rate increase was 5.2 (CI=[0.3, 10.5]; BF=47.8; P(S>C)=0.98) greater in females compared to males.

3.5 Cortisol levels

Compared to the no-stress control induction, social stress induction will increase cortisol levels (after the stress tasks).

 my_data.cortisol = 
  my_data.cortisol %>% 
  .[, concentr := mean(result), by = .(Sample.ID)] %>% 
  .[, cv.concentr := sd(result)/concentr, by = .(Sample.ID)] %>% 
  .[, .(participant, Drug, Stress, stage, plate, concentr, gender, session, time)] %>% 
  unique() %>% 
  .[, l.concentr := log(concentr)] %>% 
  .[, Sex := ifelse(gender == "f","women","men")] %>%
  .[, induced.state := ifelse(Stress == "stress" & stage > 3, "stress", "control")] %>% 
  .[, l.bl := mean(l.concentr[stage %in% 2:3]), by = .(participant, Drug, Stress)] %>% 
  .[, l.bl.m := mean(l.concentr[stage %in% 2:3], na.rm = TRUE)] %>% 
  .[, l.concentr.blc := l.concentr-l.bl+l.bl.m] %>% 
  .[, bl := mean(concentr[stage %in% 2:3]), by = .(participant, Drug, Stress)] %>% 
  .[, bl.m := mean(concentr[stage %in% 2:3], na.rm = TRUE)] %>% 
  .[, concentr.blc := concentr-bl+bl.m] %>% 
  setnames("Stress","state.condition") 

my_data.cortisol %>% 
  .[, .(N = length(unique(participant))), by = .(Drug,state.condition,stage)] %>%
  dcast(Drug+state.condition~stage, value.var = "N") %>% 
  kable(caption = "Number of participants with cortisol data by drug, state.condition and stage.") %>% kable_styling(full_width = FALSE) %>% 
  kable_classic()
Table 3.11: Number of participants with cortisol data by drug, state.condition and stage.
Drug state.condition 2 3 4 6 8 10
oxycodone control 56 56 57 58 58 58
oxycodone stress 57 60 57 56 57 58
placebo control 60 60 59 60 58 59
placebo stress 60 60 56 56 58 57
my_data.cortisol %>% 
  .[, response := concentr] %>% 
  plot_stages(my_stages = 1:10, ymin.min = NA) + 
  ylab("Cortisol concentration (mg/dl)")
Cortisol levels over time. Lines are mean responses.

Figure 3.21: Cortisol levels over time. Lines are mean responses.

my_data.cortisol %>% 
  .[, response := concentr.blc] %>% 
  plot_stages(my_stages = 1:10, ymin.min = NA) + 
  ylab("Cortisol concentration (mg/dl)")
Baseline corrected cortisol levels over time. Lines are mean responses.

Figure 3.22: Baseline corrected cortisol levels over time. Lines are mean responses.

3.5.1 By drug

my_data.cortisol %>% 
  plot_stages(my_stages = 1:10, by.Drug = TRUE, fill = "Drug") + 
  theme_classic() + 
  ylab("Plasma cortisol concentration (mg/dl)") 
Cortisol levels over time by drug. Lines are mean responses.

Figure 3.23: Cortisol levels over time by drug. Lines are mean responses.

From here on we use cortisol data collected until and including the second reminder for plotting. We do not use baseline corrected data for the statistical, because this leads to more missing data where participants do not have samples from the baseline period. In the statistical model baseline correction is achieved by adding a random intercept per participant.

my_data.cortisol[, response := concentr] 
the_data = 
  make_the_data(
    my_data.cortisol[stage < 10], n.post.stages = 10)[, outcome := "cortisol"][, session := factor(session)]
the_data %>% 
  plot_pre_post(by.Sex = FALSE, phase = "drug:induction") + 
  ylim(-1, 30) + 
  ylab("Difference plasma ortisol concentration (mg/dl)")
Change in cortisol levels from pre to post state induction. Each line represents an individual. Stress is reflected in the slope of the lines.

Figure 3.24: Change in cortisol levels from pre to post state induction. Each line represents an individual. Stress is reflected in the slope of the lines.

the_data %>% plot_pre_post(by.Sex = TRUE) +  ylim(-1, 30) + 
  ylab("Difference plasma cortisol concentration (mg/dl)") + scale_col_stress
Change in cortisol levels from pre to post state induction by sex. Each line represents an individual. Stress is reflected in the slope of the lines.

Figure 3.25: Change in cortisol levels from pre to post state induction by sex. Each line represents an individual. Stress is reflected in the slope of the lines.

p = plot_stresseff_raw(the_data, by.Sex = TRUE)
p + ggtitle("stress_effect = cortisol_in_stress-corisol_in_control")
Distribution of change in plasma cortisol levels from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Figure 3.26: Distribution of change in plasma cortisol levels from pre to post state induction, stratified by drug condition (panels) and sex (colors).

attr(p,"table") %>% 
  kable(digits = 1,
        caption = "Effects of stress by drug and sex directly calculated from raw data.") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic()
Table 3.12: Effects of stress by drug and sex directly calculated from raw data.
Drug Sex N stress_effect sd se ci
oxycodone men 24 2.3 4.0 0.8 1.7
oxycodone women 31 2.0 3.4 0.6 1.3
placebo men 26 3.0 3.9 0.8 1.6
placebo women 31 2.1 3.4 0.6 1.3
tmp = 
  p$data %>% 
  .[, .(participant,Sex,Drug, stress_effect)] %>% 
  .[, Outcome := "Plasma cortisol, ng/ML"] %>% 
  .[, stress_effect := stress_effect*10] # get on the correct scale
raw_contrasts.stress = rbind(
  raw_contrasts.stress,
  tmp
)
rm(tmp)
p
Distribution of change in plasma cortisol levels from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Figure 3.27: Distribution of change in plasma cortisol levels from pre to post state induction, stratified by drug condition (panels) and sex (colors).

plot_stresseff_raw2(the_data, by.Sex = TRUE) 
Mean and CIs of change of plasma cortisol levels (mg/dl) from pre to post state induction, stratified by drug condition and sex (colors).

Figure 3.28: Mean and CIs of change of plasma cortisol levels (mg/dl) from pre to post state induction, stratified by drug condition and sex (colors).

3.5.2 Statistical analysis: Five post induction measures

To precisely estimate the effect of stress induction on cortisol levels we estimate a hierarchical linear regression model on the cortisol levels (not baseline corrected, until including the 2nd reminder), using the model

cortisol_formula.drug = 
  bf(concentr ~ 0 + Sex:state.condition:Drug:stage + s(time, by = Sex) + (1|session) + (1 | participant) + (1 | plate) )
cortisol_formula.drug = 
  bf(concentr ~ 0 + Sex:state.condition:Drug:stage + s(time, by = Sex) + (1|session) + (state.condition:Drug:stage | participant) + (1 | plate) )

This model estimates log cortisol levels for the 16 “experimental cells” that are generated by crossing the factors state induction, drug and sex and pre/post induction. In addition, the model estimates effect of time of day and random effects for participant, stage, and session).

bfn = "brmsP1/cortisol_drug_plate+.Rdata"
if (file.exists(bfn)) {
  load(bfn) 
} else {
  my_prior = c(
    prior(normal(10,2), class = "b"),
    prior(normal(0,5), class = "b", coef = "stime:Sexman_1"),
    prior(normal(0,5), class = "b", coef = "stime:Sexwoman_1"),
    prior(normal(0,5), class = "sd", coef = "Intercept", group = "participant"),
    prior(normal(0,2), class = "sd"),
    prior(normal(0,5), class = "sds"),
    prior(normal(0,3), class = "sigma")
  )
  the_data %>% 
    .[, session := factor(session)] %>% 
    .[, participant := factor(participant)]
  bfit = 
    brm(cortisol_formula.drug,
        prior = my_prior,
        family = student(),
        data = the_data,
        iter = 2000,
        threads = threading(2, grainsize = ceiling(926/2)),
        control = list(adapt_delta = .8),
        backend = "cmdstanr")
  save(bfit, file = bfn)
}
Click to show sample diagnostics
Table 3.13: Summary of model paramters
parameter mean sd q5 q95 rhat ess_bulk ess_tail
b_Sexman:state.conditioncontrol:Drugoxycodone:stagepre.induction 9.6 0.67 8.5 10.7 1.001 1416 2434
b_Sexwoman:state.conditioncontrol:Drugoxycodone:stagepre.induction 9.0 0.64 7.9 10.0 1.005 1187 1997
b_Sexman:state.conditionstress:Drugoxycodone:stagepre.induction 9.4 0.76 8.2 10.7 1.002 1428 2556
b_Sexwoman:state.conditionstress:Drugoxycodone:stagepre.induction 8.7 0.73 7.5 9.9 1.004 1416 2352
b_Sexman:state.conditioncontrol:Drugplacebo:stagepre.induction 10.4 0.81 9.1 11.8 1.000 1659 2562
b_Sexwoman:state.conditioncontrol:Drugplacebo:stagepre.induction 9.5 0.81 8.2 10.9 1.004 1518 2044
b_Sexman:state.conditionstress:Drugplacebo:stagepre.induction 10.4 0.81 9.1 11.7 1.000 1626 2337
b_Sexwoman:state.conditionstress:Drugplacebo:stagepre.induction 11.2 0.78 9.9 12.5 1.004 1441 2346
b_Sexman:state.conditioncontrol:Drugoxycodone:stagepost.induction 8.4 0.60 7.4 9.4 1.001 1342 2170
b_Sexwoman:state.conditioncontrol:Drugoxycodone:stagepost.induction 8.1 0.58 7.1 9.0 1.007 1107 1826
b_Sexman:state.conditionstress:Drugoxycodone:stagepost.induction 10.7 0.76 9.4 11.9 1.003 1552 2342
b_Sexwoman:state.conditionstress:Drugoxycodone:stagepost.induction 9.7 0.73 8.6 10.9 1.003 1389 2353
b_Sexman:state.conditioncontrol:Drugplacebo:stagepost.induction 9.6 0.67 8.5 10.7 1.001 1705 2617
b_Sexwoman:state.conditioncontrol:Drugplacebo:stagepost.induction 8.4 0.64 7.4 9.5 1.003 1564 2223
b_Sexman:state.conditionstress:Drugplacebo:stagepost.induction 12.2 0.90 10.7 13.7 1.000 1865 1959
b_Sexwoman:state.conditionstress:Drugplacebo:stagepost.induction 12.1 0.87 10.6 13.5 1.003 1636 2724
bs_stime:Sexman_1 -3.9 5.29 -12.6 5.0 1.000 2995 2331
bs_stime:Sexwoman_1 -4.0 4.63 -11.4 3.7 1.000 5263 2849
sd_participant__Intercept 4.4 0.48 3.7 5.3 1.001 1677 2213
sd_participant__state.conditioncontrol:Drugoxycodone:stagepre.induction 0.7 0.48 0.1 1.6 1.006 301 579
sd_participant__state.conditionstress:Drugoxycodone:stagepre.induction 2.1 0.49 1.2 2.9 1.003 759 843
sd_participant__state.conditioncontrol:Drugplacebo:stagepre.induction 3.5 0.50 2.7 4.4 1.001 904 1609
sd_participant__state.conditionstress:Drugplacebo:stagepre.induction 2.8 0.55 2.0 3.7 1.002 788 1826
sd_participant__state.conditioncontrol:Drugoxycodone:stagepost.induction 0.6 0.42 0.0 1.4 1.003 299 594
sd_participant__state.conditionstress:Drugoxycodone:stagepost.induction 2.8 0.40 2.2 3.5 1.001 1045 1535
sd_participant__state.conditioncontrol:Drugplacebo:stagepost.induction 3.0 0.38 2.4 3.6 1.001 549 740
sd_participant__state.conditionstress:Drugplacebo:stagepost.induction 4.1 0.44 3.4 4.8 1.004 1094 1297
sd_plate__Intercept 1.3 0.22 0.9 1.6 1.001 2110 2790
sd_session__Intercept 0.5 0.42 0.0 1.3 1.000 1726 2263
sigma 1.0 0.09 0.9 1.2 1.003 728 1947
model: concentr ~ 0 + Sex:state.condition:Drug:stage + s(time, by = Sex) + (1 | session) + (state.condition:Drug:stage | participant) + (1 | plate)
b_ = fix effect parameters
sd_ standard deviation of hierarchical parameters
ess = effective sample size


pdata = 
  conditional_effects(bfit,"time:Sex")$time %>% 
  data.table() %>% 
  .[, .(Sex,time,estimate__,lower__,upper__)] %>% 
  .[, h := time*24] %>% 
  .[, Sex := gsub("man","men",Sex)]
pdata %>% 
  ggplot(aes(x = h, y = estimate__, fill = Sex, color = Sex)) + 
  geom_ribbon(aes(ymin = lower__, ymax = upper__), alpha = .5, color = NA) + 
  geom_line(size = 1) + 
  ylab("Plasma cortisol level") + 
  xlab("Time") + 
  scale_col_sex +
  scale_x_continuous(breaks = seq(10,20,2),
                     labels = paste0(seq(10,20,2),":00"))
estimated effect of time

Figure 3.29: estimated effect of time

The following figure shows a posterior predictive check, which confirms that the model captures drug and stress effects well.

draws = 
  bfit$fit %>% 
  as_draws() %>% 
  subset_draws("b_Sex", regex = TRUE) %>% 
  as_draws_matrix() %>% 
  data.table() %>% 
  .[, iter := 1:.N] %>% 
  melt(id.vars = "iter", variable.name = "cond", value.name = "est") %>% 
  .[grepl("pre|post",cond)] %>% 
  .[, cond := gsub("b_Sex|Drug|state.condition|stage","",cond)] %>% 
  .[, c("Sex","state.condition","Drug","stage") := tstrsplit(cond,":")] %>% 
  .[, cond := NULL] %>% 
  .[, Sex := gsub("man","men",Sex)]


obs = the_data[, .(est = mean(concentr)), by = c("Sex","state.condition","Drug","stage")]


obs %>% 
  ggplot(aes(x = stage, y=est, color = state.condition)) + 
  geom_point(aes(shape = state.condition),position = position_dodge(.2)) +
  geom_line(aes(lty = state.condition,group = factor(state.condition):factor(Drug)),
            position = position_dodge(.2)) +
  facet_grid(Drug~Sex) + 
  stat_interval(position = position_dodge(.2), data = draws, alpha = .25) +
  facet_grid(Drug~Sex) +
  scale_col_stress + 
  ylab("Cortisol level")
Observed and estimated cortisol levels pre and post state induction. Points are means of observed data and confidence bands show the 95% credible intervals of the estimated cortisol levels, after adjustment for session and time effects.

Figure 3.30: Observed and estimated cortisol levels pre and post state induction. Points are means of observed data and confidence bands show the 95% credible intervals of the estimated cortisol levels, after adjustment for session and time effects.

We use the estimated cortisol levels in the 16 cells (2 drug conditions * 2 state induction conditions * 2 sexes * 2 times) to calculate the contrasts of interest.

  1. We calculate the effect of induction as the post-pre differences in the eight cells
  2. We calculate the stress effect by subtracting the post-pre difference in the control condition from the post-pre in the stress condition, separately for Drug condition and Sex
  3. From these 4 stress effects, we can calculate average stress effects and comparisons between drug conditions or sexes
draws = 
  bfit$fit %>% 
  as_draws() %>% 
  subset_draws("b_Sex", regex = TRUE) %>% 
  as_draws_matrix() %>% 
  data.table()

setnames(draws,
         names(draws),
         gsub("b_Sex|Drug|state.condition|stage|.induction","",names(draws)))
setnames(draws,
         names(draws),
         gsub(":",".",names(draws)))
setnames(draws,names(draws),gsub("man","men", names(draws)))

make_contrasts = function(draws) {
  # put onto scale used for paper
  draws = draws * 10
  tmp = draws %>% 
    .[, `:=`(
      men.stress.oxycodone = men.stress.oxycodone.post-men.stress.oxycodone.pre,
      men.stress.placebo = men.stress.placebo.post-men.stress.placebo.pre,
      men.control.oxycodone = men.control.oxycodone.post-men.control.oxycodone.pre,
      men.control.placebo = men.control.placebo.post-men.control.placebo.pre,
      women.stress.oxycodone = women.stress.oxycodone.post-women.stress.oxycodone.pre,
      women.stress.placebo = women.stress.placebo.post-women.stress.placebo.pre,
      women.control.oxycodone = women.control.oxycodone.post-women.control.oxycodone.pre,
      women.control.placebo = women.control.placebo.post-women.control.placebo.pre
    )] %>% 
    .[, .(
      men.stress_effect.oxycodone = men.stress.oxycodone - men.control.oxycodone,
      men.stress_effect.placebo = men.stress.placebo - men.control.placebo,
      women.stress_effect.oxycodone = women.stress.oxycodone - women.control.oxycodone,
      women.stress_effect.placebo = women.stress.placebo - women.control.placebo,
      men.drug_effect.stress = men.stress.oxycodone - men.stress.placebo,
      men.drug_effect.control = men.control.oxycodone - men.control.placebo,
      women.drug_effect.stress = women.stress.oxycodone - women.stress.placebo,
      women.drug_effect.control = women.control.oxycodone - women.control.placebo
    )] %>% 
    .[,`:=`(men.drug_x_stress = men.stress_effect.oxycodone - men.stress_effect.placebo,
            women.drug_x_stress = women.stress_effect.oxycodone - women.stress_effect.placebo
    )] %>% 
    .[, iter := 1:.N] %>% 
    melt(id.vars = "iter", value.name = "est", variable.name = "contr") %>% 
    .[, Sex := ifelse(grepl("women",contr),"women","men")] %>% 
    .[, contr := gsub("women.|men.","",contr)] %>% 
    dcast(iter + contr ~ Sex, value.var = "est") %>% 
    .[, `:=`(avg = (men+women)/2, `women-men` = women-men)] %>% 
    melt(id.var = c("iter","contr"), value.name = "est", variable.name = "Sex") %>%
    dcast(iter + Sex ~ contr) %>% 
    .[, stress_effect := (stress_effect.oxycodone+stress_effect.placebo)/2] %>% 
    .[, drug_effect := (drug_effect.stress+drug_effect.control)/2]
  return(tmp)
}

contrast.samples = make_contrasts(draws)

Here are main and interaction effects of state induction and drug:

contrast.samples %>% 
  .[, .(iter,Sex,drug_effect,stress_effect,drug_x_stress)] %>% 
  melt(id.var = c("iter","Sex"), variable.name = "contrast", value.name = "est") %>% 
  .[Sex == "avg"] %>% 
  ggplot(aes(x = est)) +
  facet_wrap(~contrast, nrow = 1) +
  geom_vline(xintercept = 0) +
  stat_halfeye() + 
  coord_flip() + 
  xlab("change of cortisol level (ng/ML)") +
  ylab("Posterior density")
Stress and drugs effects on log cortisol levels

Figure 3.31: Stress and drugs effects on log cortisol levels

write.cortisol.stat = function(var) {
  contrast.samples[Sex == "avg",get(var)] %>% 
    print_effect(digits = 1, BF.thresh = 0)
}

The estimated effect of stress from the state induction on log cortisol levels in ng/ML is over all 22.9 (CI=[14.4, 31.7]; BF>10^5; P(S>C)>0.999) (cortisol is higher under stress), the effect of oxycodone is -1.7 (CI=[-9.5, 5.7]; BF=0.5; P(S>C)=0.32) (cortisol is not credibly lower under oxycodone), and the interaction effect is -1.1 (CI=[-15.3, 13.2]; BF=0.8; P(S>C)=0.44).

The next figure shows that there are no large Sex difference in the cortisol effects:

contrast.samples %>% 
  .[, .(iter,Sex,drug_effect,stress_effect,drug_x_stress)] %>% 
  melt(id.var = c("iter","Sex"), variable.name = "contrast", value.name = "est") %>% 
  .[Sex %in% c("men","women")] %>% 
  ggplot(aes(x = est, fill = Sex)) +
  facet_wrap(~contrast, nrow = 1) +
  geom_vline(xintercept = 0) +
  stat_halfeye(alpha = .5) + 
  coord_flip() + 
  xlab("change of cortisol level (ng/ML)") + 
  scale_col_sex +
  ylab("Posterior density")
Stress and drugs effects on log cortisol levels, stratified by sex.

Figure 3.32: Stress and drugs effects on log cortisol levels, stratified by sex.

Here is the table with the relevant statistics.

contrast.samples %>% 
  .[Sex == "women-men",
    .(iter,Sex,drug_effect,stress_effect,drug_x_stress)] %>% 
  melt(id.var = c("iter","Sex"), variable.name = "contrast", value.name = "sex.diff") %>% 
  .[, .(estimate = mean(sex.diff),
        CI_low = quantile(sex.diff, probs = .025),
        CI_high = quantile(sex.diff, probs = .975),
        `P > 0` = mean(sex.diff > 0)),
    by = contrast] %>% 
  kable(digits = 2,
        caption = "Sex difference in cortisol analyses. Numbers are sex difference for the contrasts in ng/ML.") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic_2() 
Table 3.14: Sex difference in cortisol analyses. Numbers are sex difference for the contrasts in ng/ML.
contrast estimate CI_low CI_high P > 0
drug_effect 6.15 -8.76 20.74 0.79
stress_effect -6.84 -24.21 10.27 0.22
drug_x_stress 1.80 -26.69 29.97 0.55
# update table with stress effects
tmp = 
  contrast.samples[, .(Sex,stress_effect)] %>% 
  .[, .(get_mci(stress_effect, get.P = ifelse(grepl("-",Sex),TRUE,FALSE))), by = .(Sex)] %>% 
  dcast(.~Sex) %>% 
  setnames(old.nms,new.nms) %>% 
  .[,new.nms,with = FALSE]
summary_table.stress = rbind(
  summary_table.stress,
  cbind(Outcome = "Plasma cortisol, ng/ML",tmp)
)
rm(tmp)

3.6 Summary of stress effects

summary_table.stress %>% 
  kable(caption = "Effects of Stress") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic() %>% 
  add_footnote("Numbers are mean [Credible interval] and posterior probability that an effect is larger than 0. The latter is only given for sex differences.")
Table 3.15: Effects of Stress
Outcome men women total Sex difference
Stressed (VAS) 29 [24, 34]; 48 [42, 53]; 38 [34, 42]; 18 [11, 26];
Negative Mood 13 [11, 16]; 26 [23, 28]; 20 [18, 21]; 12 [9, 16];
Positive Mood -20 [-23, -17]; -35 [-38, -32]; -27 [-29, -25]; -15 [-19, -11];
Heart rate, BPM 10 [6, 13]; 15 [12, 19]; 12 [10, 15]; 5 [0, 11]; 0.9795
Plasma cortisol, ng/ML 26 [14, 39]; 19 [7, 31]; 23 [14, 32]; -7 [-24, 10]; 0.21875
a Numbers are mean [Credible interval] and posterior probability that an effect is larger than 0. The latter is only given for sex differences.
summary_stats.stress = 
  summary_table.stress %>% 
  .[, .(Outcome, men, women)] %>% 
  melt(id.var = "Outcome", variable.name = "Sex") %>% 
  .[, value := gsub("\\]|;","",value)] %>% 
  .[, c("stress_effect","lower","upper") := tstrsplit(value,", | \\[")] %>% 
  .[, value := NULL] %>% 
  .[, `:=`(stress_effect = as.numeric(stress_effect), lower = as.numeric(lower), upper = as.numeric(upper))] %>% 
  .[Outcome %in% unique(raw_contrasts.stress$Outcome)]
raw_contrasts.stress %>% 
  .[, Sex := gsub("a","e",Sex)] %>% 
  .[, l := fmean(stress_effect)-3*fsd(stress_effect), by = .(Outcome,Sex)] %>% 
  .[, u := fmean(stress_effect)+3*fsd(stress_effect), by = .(Outcome,Sex)] %>% 
  .[, plot := stress_effect < u & stress_effect > l] %>% 
  .[, Drug := stringr::str_to_title(Drug)]

summary_stats.stress[, Outcome := paste("\u0394 ", Outcome)]
raw_contrasts.stress[, Outcome := paste("\u0394 ", Outcome)]
breaks_stress_fun <- function(x) {
  s.count <<- s.count + 1L
  switch(
    s.count,
    seq(0,60,20), c(0),
    seq(-50,150,50),c(0),
    seq(0,100,25), c(0),
    seq(0,60,20), c(0),
    seq(-50,150,50),c(0),
    seq(0,100,25)
  )
}
s.count = 0
summary_plot.stress = 
  summary_stats.stress %>% 
  ggplot(aes(x = Sex, y = stress_effect)) + 
  geom_bar(stat = "identity", aes(fill = Sex)) + 
  geom_hline(yintercept = 0) +
  geom_errorbar(aes(ymin = lower, ymax = upper), 
                linewidth = 1, width = 1/5) +
  facet_grid(Outcome~., scale = "free_y",switch = "y") + 
  scale_y_continuous(breaks = breaks_stress_fun) +
  scale_shape_manual(values = c(3,4)) +
  geom_quasirandom(data = raw_contrasts.stress[stress_effect < u & stress_effect > l],
                   aes(shape = Drug), color = adjustcolor("black",alpha = .5)) +
  guides(fill = guide_legend(nrow = 1, title.position = "left"), 
#         shape = guide_legend(nrow = 1, title.position = "left"),
         shape = "none") + 
  theme(legend.position = "top", axis.line.y = element_line(), 
        panel.grid.major.y = element_blank(),panel.spacing = unit(1, "lines"),
        strip.background = element_blank(), legend.key.size = unit(1,"line"),
        strip.placement = "outside", legend.title=element_text(size=12)) + 
  scale_col_sex2 + ylab("") + xlab("") + gg_no_x_axis

3.7 Reminders

3.7.1 Embarrassed (1st reminder)

We repeat the same analysis as for stress for the first reminder.

my_data.Q.emb = select_data(c("embarassed"))
my_stages = 5:6
my_data.Q.emb %>% plot_stages(my_stages = my_stages)
Embarrassed ratings over time. Lines are mean responses.

Figure 3.33: Embarrassed ratings over time. Lines are mean responses.

the_data = make_the_data(my_data.Q.emb, phase = "reminder1")[, outcome := "emb"]
attr(the_data,"N")
Table 3.16: Number of particicipants per condition by condition and stage
Drug state.condition pre.reminder1 post.reminder1
oxycodone control 62 62
oxycodone stress 63 63
placebo control 63 63
placebo stress 63 63
the_data %>% plot_pre_post("reminder") + scale_col_stress
Change in embarrassed ratings from pre to post state induction. Each line represents an individual. The stress induction effect is reflected in the slope of the lines.

Figure 3.34: Change in embarrassed ratings from pre to post state induction. Each line represents an individual. The stress induction effect is reflected in the slope of the lines.

We can further summarize this information by showing the distribution of the slopes in the previous figure:

plot_stresseff_raw(the_data)
Distribution on change in embarrassed ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Figure 3.35: Distribution on change in embarrassed ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

plot_stresseff_raw2(the_data)
Mean and CIs of change in embarrassed ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Figure 3.36: Mean and CIs of change in embarrassed ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Here we run the Bayesian ordinal logistic regression:

fn = "brmsP1/Reminder/emb1_acat.Rdata"
stress_formula = response.c ~ stage:induced.state:Sex + session + (1|participant) 
the_data = the_data %>% 
  coarsen(var = "response", bins = 20, range = c(0,100))
if (file.exists(fn)) {
  load(fn)
} else {
  fit = fit_ordinal_model(the_data, stress_formula, fn, family = acat())
}
#pp_check(fit, ndraws = 2000, type = "bars")
Click to show sample diagnostics
Table 3.17: Summary of model paramters
parameter mean sd q5 q95 rhat ess_bulk ess_tail
Threshold[1] 1.6 0.71 0.4 2.8 1.003 1066 1637
Threshold[2] 0.7 0.72 -0.5 1.9 1.005 1102 1720
Threshold[3] 0.4 0.74 -0.8 1.6 1.002 1169 1770
Threshold[4] 0.7 0.77 -0.6 2.0 1.002 1198 1912
Threshold[5] 0.8 0.79 -0.5 2.0 1.005 1217 1977
Threshold[6] 0.9 0.82 -0.5 2.3 1.001 1317 2212
Threshold[7] 1.0 0.90 -0.4 2.6 1.002 1518 2177
Threshold[8] 0.8 0.95 -0.7 2.4 1.003 1556 2354
Threshold[9] -0.3 0.91 -1.8 1.1 1.001 1601 2542
Threshold[10] 2.0 0.93 0.5 3.5 1.000 1615 2377
Threshold[11] 0.7 1.05 -1.0 2.4 1.001 1636 2216
Threshold[12] 0.0 0.99 -1.7 1.6 1.001 1951 2550
Threshold[13] 0.8 0.90 -0.7 2.3 1.002 1530 2381
Threshold[14] 0.5 0.87 -0.9 2.0 1.004 1442 2194
Threshold[15] 3.1 1.14 1.3 5.1 1.001 2045 2822
Threshold[16] 0.2 1.25 -2.0 2.2 1.001 2040 2630
Threshold[17] 0.5 1.08 -1.3 2.2 1.002 2027 2547
Threshold[18] 0.4 1.00 -1.3 2.0 1.000 1811 2436
b_session2 0.1 0.07 0.0 0.2 1.001 2641 3169
b_session3 0.1 0.05 0.0 0.2 1.000 4649 3033
b_session4 0.1 0.07 0.0 0.2 1.000 3124 3268
b_stagepre.reminder1:induced.statecontrol:Sexman -0.2 0.70 -1.4 0.9 1.004 1034 1569
b_stagepost.reminder1:induced.statecontrol:Sexman -0.4 0.70 -1.5 0.8 1.003 1038 1613
b_stagepre.reminder1:induced.statestress:Sexman 0.0 0.69 -1.2 1.1 1.004 1020 1515
b_stagepost.reminder1:induced.statestress:Sexman 0.2 0.69 -0.9 1.4 1.004 1025 1660
b_stagepre.reminder1:induced.statecontrol:Sexwoman -0.1 0.70 -1.2 1.1 1.004 1028 1736
b_stagepost.reminder1:induced.statecontrol:Sexwoman -0.5 0.70 -1.6 0.7 1.004 1029 1527
b_stagepre.reminder1:induced.statestress:Sexwoman 0.3 0.69 -0.8 1.4 1.005 1017 1652
b_stagepost.reminder1:induced.statestress:Sexwoman 0.6 0.69 -0.5 1.7 1.005 1016 1639
sd_participant__Intercept 0.3 0.06 0.2 0.5 1.002 1196 2300
model: response.c ~ stage:induced.state:Sex + session + (1 | participant)
b_ = fix effect parameters
sd_ standard deviation of hierarchical parameters
ess = effective sample size


We calculate the contrasts manually from posterior predictions:

tmp = plot_prepost_contrast(fit, the_data,"emb1", my_contrast = "stress")
tmp[[1]]
Posterior distributions of estimated first stress reminder effects on self reported embarrassment on the 0-100 VAS scale in women, men, and sex difference (women-men). Density plots show posterior distributions. Dots and lines identify means, and 50 & 95% credible intervals.

Figure 3.37: Posterior distributions of estimated first stress reminder effects on self reported embarrassment on the 0-100 VAS scale in women, men, and sex difference (women-men). Density plots show posterior distributions. Dots and lines identify means, and 50 & 95% credible intervals.

tmp$table
Table 3.18: Modeled post minus pre state induction on emb1 ratings by condition (control, stress) and the stress effect by sex.
Sex control stress stress_effect
men -1 [-4, 1]; 8 [3, 12]; 9 [4, 14];
women -3 [-6, -1]; 20 [13, 26]; 23 [16, 29];
total -2 [-4, -0]; 14 [10, 18]; 16 [11, 20];
Sex difference -2 [-5, 2]; 12 [4, 19]; 14 [5, 22];
post = tmp[[2]]

The estimated effect of stress from the first state reminder is NaN (CI=[NA, NA]; BFNA; P(S>C)NA). The effect is 23 (CI=[16, 29]; BF>10^5; P(S>C)>0.999) in females, 9 (CI=[4, 14]; BF=1999; P(S>C)>0.999) in males and the sex difference is NaN (CI=[NA, NA]; BFNA; P(S>C)NA).

3.7.2 Embarrassed

We repeat the same analysis as for stress and also take into account that there are multiple reminders.

my_data.Q.emb = select_data(c("embarassed"))
my_stages = 5:8
my_data.Q.emb %>% plot_stages(my_stages = my_stages)
Embarrassed ratings over time. Lines are mean responses.

Figure 3.38: Embarrassed ratings over time. Lines are mean responses.

the_data = make_the_data(my_data.Q.emb, phase = "reminder")[, outcome := "emb"]
attr(the_data,"N")
Table 3.19: Number of particicipants per condition by condition and stage
Drug state.condition pre.reminder post.reminder
oxycodone control 62 62
oxycodone stress 63 63
placebo control 63 63
placebo stress 63 63
the_data %>% plot_pre_post("reminder") + scale_col_stress
Change in embarrassed ratings from pre to post state induction. Each line represents an individual. The stress induction effect is reflected in the slope of the lines.

Figure 3.39: Change in embarrassed ratings from pre to post state induction. Each line represents an individual. The stress induction effect is reflected in the slope of the lines.

We can further summarize this information by showing the distribution of the slopes in the previous figure:

plot_stresseff_raw(the_data)
Distribution on change in embarrassed ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Figure 3.40: Distribution on change in embarrassed ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

plot_stresseff_raw2(the_data)
Mean and CIs of change in embarrassed ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Figure 3.41: Mean and CIs of change in embarrassed ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Here we run the Bayesian ordinal logistic regression:

fn = "brmsP1/Reminder/emb_acat.Rdata"
stress_formula = response.c ~ induced.state + Sex + induced.state:Sex + session + (1|participant) + (1|stage)
the_data = the_data %>% 
  coarsen(var = "response", bins = 20, range = c(0,100))
if (file.exists(fn)) {
  load(fn)
} else {
  fit = fit_ordinal_model(the_data, stress_formula, fn, family = acat())
}
#pp_check(fit, ndraws = 2000, type = "bars")
Click to show sample diagnostics
Table 3.20: Summary of model paramters
parameter mean sd q5 q95 rhat ess_bulk ess_tail
Threshold[1] 1.8 0.28 1.2 2.2 1.015 427 349
Threshold[2] 0.8 0.30 0.3 1.3 1.013 458 340
Threshold[3] 0.7 0.33 0.2 1.2 1.009 519 400
Threshold[4] 0.5 0.34 -0.1 1.1 1.008 555 511
Threshold[5] 0.8 0.37 0.2 1.4 1.005 754 778
Threshold[6] 1.1 0.41 0.4 1.8 1.007 665 469
Threshold[7] 1.0 0.50 0.2 1.9 1.001 1109 1415
Threshold[8] 0.5 0.52 -0.4 1.3 1.004 1148 1678
Threshold[9] 0.4 0.47 -0.4 1.1 1.004 987 1474
Threshold[10] 0.9 0.48 0.1 1.6 1.003 789 497
Threshold[11] 0.5 0.47 -0.3 1.3 1.003 954 1305
Threshold[12] 0.7 0.46 0.0 1.5 1.003 1112 1955
Threshold[13] 1.4 0.54 0.5 2.3 1.003 1357 1638
Threshold[14] 0.3 0.55 -0.6 1.2 1.002 1334 1980
Threshold[15] 1.9 0.63 0.9 3.0 1.003 1459 1893
Threshold[16] 0.8 0.73 -0.4 2.0 1.003 884 531
Threshold[17] 0.6 0.71 -0.6 1.8 1.001 1840 1558
Threshold[18] 1.4 0.75 0.2 2.7 1.004 1152 975
Threshold[19] -0.6 0.68 -1.8 0.4 1.000 1708 1889
b_induced.statestress 0.4 0.06 0.3 0.5 1.003 1550 951
b_Sexwoman 0.2 0.11 0.0 0.3 1.012 696 1291
b_session2 0.0 0.04 -0.1 0.1 1.003 1640 2245
b_session3 0.0 0.03 0.0 0.1 1.003 3359 2516
b_session4 -0.1 0.04 -0.2 0.0 1.002 1559 2418
b_induced.statestress:Sexwoman 0.1 0.07 0.0 0.2 1.002 2012 1613
sd_participant__Intercept 0.3 0.05 0.2 0.4 1.008 675 1675
sd_stage__Intercept 0.4 0.35 0.1 1.1 1.017 386 1724
model: response.c ~ induced.state + Sex + induced.state:Sex + session + (1 | participant) + (1 | stage)
b_ = fix effect parameters
sd_ standard deviation of hierarchical parameters
ess = effective sample size


We calculate the contrasts manually from posterior predictions:

tmp = plot_prepost_contrast(fit, the_data,"emb", my_contrast = "stress")
tmp[[1]]
Posterior distributions of estimated stress reminders effects on self reported embarrassment on the 0-100 VAS scale in women, men, and sex difference (women-men). Density plots show posterior distributions. Dots and lines identify means, and 50 & 95% credible intervals.

Figure 3.42: Posterior distributions of estimated stress reminders effects on self reported embarrassment on the 0-100 VAS scale in women, men, and sex difference (women-men). Density plots show posterior distributions. Dots and lines identify means, and 50 & 95% credible intervals.

tmp$table
Table 3.21: Modeled post minus pre state induction on emb ratings by condition (control, stress) and the stress effect by sex.
Sex control stress stress_effect
men 1 [1, 2]; 7 [5, 9]; 5 [4, 8];
women 2 [1, 3]; 14 [10, 18]; 12 [9, 16];
total 2 [1, 2]; 11 [8, 14]; 9 [6, 12];
Sex difference 0 [-0, 1]; 7 [5, 10]; 7 [4, 10];
post = tmp[[2]]

The estimated effect of stress from the state reminder on feeling embarrased is NaN (CI=[NA, NA]; BFNA; P(S>C)NA). The effect is 12 (CI=[9, 16]; BF>10^5; P(S>C)>0.999) in females, 5 (CI=[4, 8]; BF>10^5; P(S>C)>0.999) in males and the sex difference is NaN (CI=[NA, NA]; BFNA; P(S>C)NA).

3.7.3 Selfconscious

We repeat the same analysis as for stress and also take into account that there are multiple reminders.

my_data.Q.selfc = select_data(c("selfconscious"))
my_stages = 5:8
my_data.Q.selfc %>% plot_stages(my_stages=my_stages)
Embarrassed ratings over time. Lines are mean responses.

Figure 3.43: Embarrassed ratings over time. Lines are mean responses.

the_data = make_the_data(my_data.Q.selfc, phase = "reminder")[, outcome := "selfc"]
attr(the_data,"N")
Table 3.22: Number of particicipants per condition by condition and stage
Drug state.condition pre.reminder post.reminder
oxycodone control 62 62
oxycodone stress 63 63
placebo control 63 63
placebo stress 63 63
the_data %>% plot_pre_post("reminder") + scale_col_stress
Change in selfconscious ratings from pre to post state induction. Each line represents an individual. The stress induction effect is reflected in the slope of the lines.

Figure 3.44: Change in selfconscious ratings from pre to post state induction. Each line represents an individual. The stress induction effect is reflected in the slope of the lines.

We can further summarize this information by showing the distribution of the slopes in the previous figure:

plot_stresseff_raw(the_data)
Distribution on change in embarrassed ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Figure 3.45: Distribution on change in embarrassed ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

plot_stresseff_raw2(the_data)
Mean and CIs of change in selfconcious ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Figure 3.46: Mean and CIs of change in selfconcious ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Here we run the Bayesian ordinal logistic regression:

fn = "brmsP1/Reminder/selfc_acat.Rdata"
stress_formula = response.c ~ induced.state + Sex + induced.state:Sex + session + (1|participant) + (1|stage)
the_data = the_data %>% 
  coarsen(var = "response", bins = 20, range = c(0,100))
if (file.exists(fn)) {
  load(fn)
} else {
  fit = fit_ordinal_model(the_data, stress_formula, fn, family = acat())
}
#pp_check(fit, ndraws = 2000, type = "bars")

We calculate the contrasts manually from posterior predictions:

tmp = plot_prepost_contrast(fit, the_data,"selfc", my_contrast = "stress")
tmp[[1]]
Stress effects in women, men, and sex difference (women-men)

Figure 3.47: Stress effects in women, men, and sex difference (women-men)

tmp$table
Table 3.23: Modeled post minus pre state induction on selfc ratings by condition (control, stress) and the stress effect by sex.
Sex control stress stress_effect
men 1 [-0, 3]; 1 [-0, 3]; -0 [-0, 0];
women 2 [-0, 4]; 2 [-0, 4]; -0 [-0, 0];
total 2 [-0, 4]; 2 [-0, 4]; -0 [-0, 0];
Sex difference 0 [-0, 1]; 0 [-0, 1]; 0 [-0, 0];
post = tmp[[2]]

The estimated effect of stress from the state reminder on feeling selfconscious is NaN (CI=[NA, NA]; BFNA; P(S>C)NA). The effect is 0 (CI=[0, 0]; BF=0.1; P(S>C)=0.06) in females, 0 (CI=[0, 0]; BF=0.1; P(S>C)=0.06) in males and the sex difference is NaN (CI=[NA, NA]; BFNA; P(S>C)NA).

3.7.4 Stressed

We repeat the same analysis as for stress and also take into account that there are multiple reminders.

my_data.Q.selfc = select_data(c("stressed"))
my_stages = 5:8
my_data.Q.selfc %>% plot_stages(my_stages=my_stages)
Embarrassed ratings over time. Lines are mean responses.

Figure 3.48: Embarrassed ratings over time. Lines are mean responses.

the_data = make_the_data(my_data.Q.selfc, phase = "reminder")[, outcome := "selfc"]
attr(the_data,"N")
Table 3.24: Number of particicipants per condition by condition and stage
Drug state.condition pre.reminder post.reminder
oxycodone control 62 62
oxycodone stress 63 63
placebo control 63 63
placebo stress 63 63
the_data %>% plot_pre_post("reminder") + scale_col_stress
Change in 'stressed' ratings from pre to post state induction. Each line represents an individual. The stress induction effect is reflected in the slope of the lines.

Figure 3.49: Change in ‘stressed’ ratings from pre to post state induction. Each line represents an individual. The stress induction effect is reflected in the slope of the lines.

We can further summarize this information by showing the distribution of the slopes in the previous figure:

plot_stresseff_raw(the_data)
Distribution on change in embarrassed ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Figure 3.50: Distribution on change in embarrassed ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

plot_stresseff_raw2(the_data)
Mean and CIs of change in selfconcious ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Figure 3.51: Mean and CIs of change in selfconcious ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Here we run the Bayesian ordinal logistic regression:

fn = "brmsP1/Reminder/stressed_acat.Rdata"
stress_formula = response.c ~ induced.state + Sex + induced.state:Sex + session + (1|participant) + (1|stage)
the_data = the_data %>% 
  coarsen(var = "response", bins = 20, range = c(0,100))
if (file.exists(fn)) {
  load(fn)
} else {
  fit = fit_ordinal_model(the_data, stress_formula, fn, family = acat())
}
#pp_check(fit, ndraws = 2000, type = "bars")

We calculate the contrasts manually from posterior predictions:

tmp = plot_prepost_contrast(fit, the_data,"stressed", my_contrast = "stress")
tmp[[1]]
Stress effects in women, men, and sex difference (women-men)

Figure 3.52: Stress effects in women, men, and sex difference (women-men)

tmp$table
Table 3.25: Modeled post minus pre state induction on stressed ratings by condition (control, stress) and the stress effect by sex.
Sex control stress stress_effect
men 1 [0, 1]; 2 [0, 4]; 1 [0, 3];
women 1 [0, 3]; 3 [0, 6]; 2 [0, 3];
total 1 [0, 2]; 2 [0, 5]; 1 [0, 3];
Sex difference 1 [0, 2]; 1 [0, 2]; 0 [-0, 1];
post = tmp[[2]]

The estimated effect of stress from the state reminder on feeling stressed is NaN (CI=[NA, NA]; BFNA; P(S>C)NA). The effect is 2 (CI=[0, 3]; BF=41.1; P(S>C)=0.98) in females, 1 (CI=[0, 3]; BF=41.1; P(S>C)=0.98) in males and the sex difference is NaN (CI=[NA, NA]; BFNA; P(S>C)NA).

3.8 Drug effects

3.8.1 High

my_data.Q.high = select_data(c("high"))
my_stages = 1:11
my_data.Q.high %>% 
  .[stage %in% my_stages] %>% 
    .[,.(m = mean(response), se = sd(response)/sqrt(length(unique(participant)))), 
      by = .(stage,Sex,Drug)] %>% 
    .[, `:=`(lower = m-2*se, upper = m+2*se)] %>% 
    ggplot(aes(x = stage, y = m, color = Drug, lty = Sex, fill = Drug)) + 
    geom_line() + 
    geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .2, color = NA) + 
    scale_x_continuous(breaks = my_stages, labels = stage.names[my_stages]) +
  scale_col_drug
Feeling high over time. Lines are mean responses.

Figure 3.53: Feeling high over time. Lines are mean responses.

the_data = make_the_data(my_data.Q.high[stage %in% 4:5],phase = "drug")[, outcome := "high"]
attr(the_data,"N")
Table 3.26: Number of particicipants per condition by condition and stage
Drug state.condition pre.drug post.drug
oxycodone control 62 62
oxycodone stress 63 63
placebo control 63 63
placebo stress 63 63
the_data %>% plot_pre_post("drug")
Change in high from drug ratings from pre to post the first drug administration. Each line represents an individual. The stress induction effect is reflected in the slope of the lines.

Figure 3.54: Change in high from drug ratings from pre to post the first drug administration. Each line represents an individual. The stress induction effect is reflected in the slope of the lines.

We can further summarize this information by showing the distribution of the slopes in the previous figure:

plot_drug_effect_raw(the_data)
Distribution on change in high from drug ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Figure 3.55: Distribution on change in high from drug ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Here we run the Bayesian ordinal logistic regression:

fn = "brmsP1/high_acat.Rdata"
the_data = the_data %>%
  .[, Drug := factor(Drug, levels = c("placebo","oxycodone"))] %>% 
  coarsen(var = "response", bins = 20, range = c(0,100))
drug_formula = response.c ~ Drug:stage:Sex + (1 | session) + (1 | participant) 
#drug_formula = response.c ~ Drug:stage:Sex + (1 | session) + (1 + Drug:stage | participant)
if (file.exists(fn)) {
  load(fn)
} else {
  fit = fit_ordinal_model(the_data, drug_formula, fn, family = acat())
}
#pp_check(fit, ndraws = 2000, type = "bars")
Click to show sample diagnostics
Table 3.27: Summary of model paramters
parameter mean sd q5 q95 rhat ess_bulk ess_tail
Threshold[1] 1.8 0.73 0.6 3.0 1.006 921 1711
Threshold[2] 0.2 0.76 -1.0 1.4 1.006 1017 1870
Threshold[3] 1.0 0.81 -0.3 2.3 1.005 1108 1948
Threshold[4] 0.5 0.88 -0.9 2.0 1.003 1255 1892
Threshold[5] -0.2 0.90 -1.6 1.3 1.003 1313 2200
Threshold[6] 0.3 0.88 -1.2 1.7 1.004 1259 2240
Threshold[7] 1.4 0.96 -0.2 3.0 1.004 1516 2189
Threshold[8] -0.1 1.03 -1.8 1.5 1.002 1578 2244
Threshold[9] -0.4 0.88 -1.9 1.0 1.003 1319 2340
Threshold[10] 1.3 0.87 -0.1 2.7 1.004 1291 2150
Threshold[11] -0.2 0.86 -1.6 1.2 1.005 1246 2295
Threshold[12] 1.0 0.84 -0.3 2.4 1.004 1276 2146
Threshold[13] 0.5 0.84 -0.9 1.9 1.004 1195 2080
Threshold[14] 1.2 0.86 -0.2 2.6 1.004 1316 2253
Threshold[15] 0.7 0.87 -0.8 2.1 1.004 1248 2289
Threshold[16] 3.1 1.16 1.3 5.0 1.002 1848 2491
Threshold[17] -1.4 1.15 -3.3 0.5 1.002 1903 2300
Threshold[18] 1.6 0.93 0.1 3.2 1.002 1526 2282
Threshold[19] -0.7 0.91 -2.3 0.7 1.005 1351 2256
b_Drugplacebo:stagepre.drug:Sexman -0.2 0.71 -1.3 1.0 1.006 872 1579
b_Drugoxycodone:stagepre.drug:Sexman -0.2 0.71 -1.4 0.9 1.006 885 1673
b_Drugplacebo:stagepost.drug:Sexman -0.2 0.71 -1.3 1.0 1.006 881 1465
b_Drugoxycodone:stagepost.drug:Sexman 0.7 0.71 -0.4 1.9 1.006 900 1637
b_Drugplacebo:stagepre.drug:Sexwoman -0.4 0.71 -1.6 0.7 1.007 892 1558
b_Drugoxycodone:stagepre.drug:Sexwoman -0.4 0.72 -1.5 0.8 1.007 895 1377
b_Drugplacebo:stagepost.drug:Sexwoman -0.2 0.71 -1.3 1.0 1.007 880 1594
b_Drugoxycodone:stagepost.drug:Sexwoman 0.8 0.71 -0.4 1.9 1.007 905 1439
sd_participant__Intercept 0.3 0.05 0.2 0.4 1.000 1591 2265
sd_session__Intercept 0.2 0.22 0.1 0.6 1.001 860 809
model: response.c ~ Drug:stage:Sex + (1 | session) + (1 | participant)
b_ = fix effect parameters
sd_ standard deviation of hierarchical parameters
ess = effective sample size


We calculate the contrasts manually from posterior predictions:

tmp = plot_prepost_contrast(fit, the_data, "high", my_contrast = "drug")
tmp[[1]]
Stress effects in women, men, and sex difference (women-men)

Figure 3.56: Stress effects in women, men, and sex difference (women-men)

tmp$table
Table 3.28: Modeled post minus pre drug administration on high ratings by condition (placebo, oxycodone) and the drug effect by sex.
Sex placebo oxycodone drug_effect
men 1 [-3, 4]; 49 [43, 55]; 49 [42, 55];
women 2 [0, 5]; 58 [53, 63]; 56 [50, 62];
total 2 [-1, 4]; 54 [50, 58]; 52 [48, 57];
Sex difference 2 [-3, 6]; 9 [1, 17]; 7 [-1, 16];
post = tmp[[2]]

The estimated effect of drug on ratings of being high is NaN (CI=[NA, NA]; BFNA; P(D>P)NA). The effect is 56 (CI=[50, 62]; BF>10^5; P(D>P)>0.999) in females, 49 (CI=[42, 55]; BF>10^5; P(D>P)>0.999) in males and the sex difference is NaN (CI=[NA, NA]; BFNA; P(D>P)NA).

3.9 Drug side effects

3.9.1 Side effect: dizzy

Raw data for ‘dizzy’

Estimated effects of Drug for ‘dizzy’

(#tab:side_effects)Modeled post minus pre drug administration on dizzy ratings by condition (placebo, oxycodone) and the drug effect by sex.
Sex placebo oxycodone drug_effect
men 2 [-1, 5]; 34 [28, 40]; 32 [26, 39];
women -1 [-6, 4]; 59 [53, 64]; 60 [52, 67];
total 0 [-2, 3]; 46 [42, 50]; 46 [41, 51];
Sex difference -3 [-9, 3]; 24 [16, 32]; 27 [18, 37];

3.9.2 Side effect: dry_mouth

Raw data for ‘dry_mouth’

Estimated effects of Drug for ‘dry_mouth’

(#tab:side_effects)Modeled post minus pre drug administration on dry_mouth ratings by condition (placebo, oxycodone) and the drug effect by sex.
Sex placebo oxycodone drug_effect
men -6 [-12, -1]; -5 [-10, -0]; 1 [-6, 8];
women -2 [-7, 4]; -1 [-6, 5]; 1 [-7, 9];
total -4 [-8, -0]; -3 [-7, 1]; 1 [-4, 6];
Sex difference 5 [-3, 12]; 4 [-3, 12]; -0 [-11, 10];

3.9.3 Side effect: blunted

Raw data for ‘blunted’

Estimated effects of Drug for ‘blunted’

(#tab:side_effects)Modeled post minus pre drug administration on blunted ratings by condition (placebo, oxycodone) and the drug effect by sex.
Sex placebo oxycodone drug_effect
men -1 [-7, 4]; 34 [27, 40]; 36 [27, 44];
women 5 [-0, 9]; 49 [43, 55]; 45 [37, 52];
total 2 [-2, 5]; 42 [37, 46]; 40 [34, 45];
Sex difference 6 [-1, 13]; 15 [6, 24]; 9 [-2, 20];

3.9.4 Side effect: nauseous

Raw data for ‘nauseous’

Estimated effects of Drug for ‘nauseous’

(#tab:side_effects)Modeled post minus pre drug administration on nauseous ratings by condition (placebo, oxycodone) and the drug effect by sex.
Sex placebo oxycodone drug_effect
men -2 [-5, 1]; 4 [1, 8]; 6 [2, 11];
women -2 [-5, 0]; 7 [3, 11]; 9 [4, 14];
total -2 [-4, -0]; 5 [3, 8]; 7 [4, 11];
Sex difference -0 [-4, 4]; 2 [-3, 8]; 3 [-4, 9];

3.9.5 Side effect: warm_face

Raw data for ‘warm_face’

Estimated effects of Drug for ‘warm_face’

(#tab:side_effects)Modeled post minus pre drug administration on warm_face ratings by condition (placebo, oxycodone) and the drug effect by sex.
Sex placebo oxycodone drug_effect
men -6 [-11, -2]; 13 [6, 20]; 19 [11, 27];
women -11 [-16, -5]; 19 [12, 26]; 30 [21, 38];
total -8 [-12, -5]; 16 [11, 21]; 24 [18, 30];
Sex difference -5 [-12, 2]; 6 [-4, 16]; 11 [-2, 22];

3.9.6 Side effect: not_myself

Raw data for ‘not_myself’

Estimated effects of Drug for ‘not_myself’

(#tab:side_effects)Modeled post minus pre drug administration on not_myself ratings by condition (placebo, oxycodone) and the drug effect by sex.
Sex placebo oxycodone drug_effect
men -2 [-7, 3]; 20 [13, 27]; 22 [14, 31];
women -5 [-11, 1]; 34 [27, 41]; 40 [30, 49];
total -4 [-8, 0]; 27 [22, 32]; 31 [25, 37];
Sex difference -3 [-11, 5]; 14 [4, 24]; 17 [4, 30];

3.10 Participants will use vigorous button presses

We display the number of button presses participants made, stratified by drug condition and target dose. Participants with a target dose larger than 50% had to press buttons to achieve a high dose and with a target dose lower than 50% had to press buttons to achieve a low dose.

dt.presses = 
  my_data %>% 
  .[cheater == FALSE] %>% 
  setkeyv(c("participant","session")) %>% 
  .[, last.trial := ifelse(time == max(time),TRUE,FALSE), by = .(session, participant)] %>% 
  .[, n.presses := p_down + p_up] %>% 
  .[, n.presses := cumsum(n.presses), by = .(session, participant)] %>% 
  .[, presses.per.s := n.presses/time] %>% 
  .[last.trial == TRUE] %>% 
  .[!is.na(satisfact)] %>% 
  .[, stage := 0] %>% 
  .[, N.oxy := sum(Drug == "oxycodone"), by = .(participant)] %>% 
  .[N.oxy > 1] %>% 
  .[, .(Sex, participant, session, Drug, 
        state.condition, target, n.presses, presses.per.s)] 


dt.presses[, `target effect` := ifelse(target > 50, "target > 50%","target < 50%")]

hist.presses = 
  dt.presses %>% 
  ggplot(aes(x = presses.per.s, fill = Drug)) + 
  geom_histogram(position = "identity", alpha = .5, bins = 15) + 
  facet_grid(Sex~`target effect`) + 
  xlab("Number of presses second") + 
  ylab("Number of participants") +
  scale_col_drug

hist.presses

dt.presses %>% 
  .[, .(mean = round(mean(presses.per.s),1), sd = round(sd(presses.per.s),1), N = length(unique(participant))),
   by = .(Drug,state.condition)] %>% 
  kable(caption = "Number of button presses per second by drug and stress condition and sex") %>% 
  kable_styling(full_width = FALSE)
Table 3.29: Number of button presses per second by drug and stress condition and sex
Drug state.condition mean sd N
placebo stress 5.6 0.5 52
placebo control 5.5 0.8 50
oxycodone stress 5.6 0.6 56
oxycodone control 5.5 0.6 56
dt.presses %>% 
  .[, .(mean = round(mean(n.presses)), sd = round(sd(n.presses))),
   by = .(Drug,state.condition,Sex)] %>% 
  kable(caption = "Number of button presses by drug and stress condition and sex") %>% 
  kable_styling(full_width = FALSE)
Table 3.30: Number of button presses by drug and stress condition and sex
Drug state.condition Sex mean sd
placebo stress women 662 64
placebo control women 651 71
oxycodone stress women 636 50
oxycodone control women 630 84
placebo control men 672 123
placebo stress men 687 69
oxycodone control men 689 51
oxycodone stress men 705 75

4 Drug wanting and liking after drug administration

The data for the second drug administration includes only participants who wanted a second dosage larger than zero.

my_data[is.na(satisfact), last.trial := ifelse(time == max(time),TRUE,FALSE), by = .(session, participant)]

my_data_beh = 
  my_data[last.trial == TRUE & is.na(satisfact)] %>% 
  .[, stage := 0] %>% 
  .[, N.oxy := sum(Drug == "oxycodone"), by = .(participant)] %>% 
  .[N.oxy > 1] %>% 
  .[, s.effect_obtained := effect_obtained-mean(effect_obtained), by = .(participant)] %>% 
  .[, s.effect_obtained := s.effect_obtained + mean(effect_obtained)] 
pdata = 
  my_data.Q %>% 
  setkeyv(key(my_data_beh)) %>% 
  .[item_name %in% c("drug_like","drug_again","drug_dislike") & stage %in% c(5,11)] %>% 
  .[my_data, effect_obtained := effect_obtained] %>% 
  .[stage == 5, effect_obtained := 100] %>% 
  .[stage == 11 & effect_obtained == 0, effect_obtained := NA] %>% 
  .[, .(response = mean(response, na.rm = TRUE), sd = sd(response, na.rm = TRUE), N = .N),
    by = .(state.condition,Drug,stage,item_name)] %>% 
  .[, se := sd/sqrt(N)] %>% 
  .[, `:=`(lower=response-2*se,
           upper=response+2*se)] %>% 
  .[, time := stage-4] %>% 
  .[time>1, time := time-4] %>% 
  .[state.condition == "control", time:=time-.05] %>% 
  .[state.condition != "control", time:=time+.05] %>% 
  .[, Stage := factor(stage.names[stage],levels = stage.names)]

pdata %>% 
  ggplot(aes(x = Stage, y = response, ymin=lower, ymax=upper,
             color = Drug, shape = state.condition), group = state.condition) +
  geom_pointrange(size = 1, position = position_dodge(.2)) +
  #geom_vline(xintercept = c(0.75,2.75), lty = 2, alpha = .5) +
  #scale_x_continuous(breaks = c(0.75,2.75), labels = c("dose 1","dose 2")) +  
  ylab("Response") + 
  guides(color=guide_legend(nrow=2,byrow=TRUE),
         shape=guide_legend(nrow=2,byrow=TRUE)) +
  scale_col_drug + 
  xlab("Time in the experiment") +
  facet_wrap(~item_name) +
  theme(legend.position = "bottom") + 
  theme(axis.line = element_line())
Drug liking, disliking, and wanting after the administration of the first and second dose.

Figure 4.1: Drug liking, disliking, and wanting after the administration of the first and second dose.

5 Main Results

5.1 Effect obtained

Prepare data

my_data[is.na(satisfact), last.trial := ifelse(time == max(time),TRUE,FALSE), by = .(session, participant)]

my_data %>% 
  .[, n.presses := p_down + p_up] %>% 
  .[, n.presses := cumsum(n.presses), by = .(session, participant)]

my_data_beh = 
  my_data[last.trial == TRUE & is.na(satisfact)] %>% 
  .[, stage := 0] %>% 
  .[, N.oxy := sum(Drug == "oxycodone"), by = .(participant)] %>% 
  .[N.oxy > 1] %>% 
  .[, s.effect_obtained := effect_obtained-mean(effect_obtained), by = .(participant)] %>% 
  .[, s.effect_obtained := s.effect_obtained + mean(effect_obtained)] 

load("data/dosage_data.Rdata")
feel_tresh = dosage_data[Drug == "placebo", quantile(drug_feel, .95)]
dosage_data = 
  dosage_data %>% 
  .[, session := as.integer(session)] %>% 
  setkeyv(.,key(my_data_beh)) %>% 
  .[my_data_beh, target := target] %>% 
  .[, drug_2_refused := 
      ifelse(is.na(drug_2_admin) & 
               drug_2_dose == 0 & 
               target > 0 & 
               drug_feel < feel_tresh, TRUE, FALSE)]

dosage.factor = 
  median(dosage_data$weight_kg/dosage_data$drug_1_dose,
       na.rm = TRUE)

mg_per_kg = round(mean(dosage_data$drug_1_dose/dosage_data$weight_kg),3)

m.format = function(x) round(mean(x,na.rm = TRUE),2)

5.1.1 Drug dosages

The first oxycodone dose was 0.043 mg per kg body weight, and an average dose of 3.12 mg. The dose for the second drug administration was determined by the result of the self administration task, whereby participant who had worked to obtain more than 100% of the effect after the first dose still received 100%. In the oxycodone condition after stress induction participants received on average 2.34 mg at the second administration (working for on average 85.02%), whereas they received after control induction on average 2.16 mg (working for on average 77.68%). Considering only those who wanted a non-zero dose, participants received on average 2.88 mg n the oxycodone condition after stress induction at the second administration (working for on average 97.4%), whereas they received after control induction on average 2.62 mg (working for on average 92.37%).

3 participants, who had worked to obtain on average 3.75% of the first effect, decided that they did not want the second drug administration. 13 participants worked to not obtain any drug at the second administration.

Drug effects split by Sex:

my_data_beh.stat = 
  my_data_beh %>% 
  .[, .(effect_obtained = mean(effect_obtained),
        se = sd(effect_obtained)/sqrt(.N),
        N = .N),
    by = .(state.condition,Drug,stage,Sex)] %>% 
  .[, `:=`(lower = effect_obtained-2*se,upper = effect_obtained+2*se)]

my_data_beh.stat %>% 
  ggplot(aes(x = state.condition, y = effect_obtained, color = Drug, group = Drug)) + 
  geom_quasirandom(data = my_data_beh, alpha = .25) + 
  geom_line(position = position_dodge(.5)) +
  geom_point(position = position_dodge(.5), size = 4) + 
  geom_linerange(aes(ymin = lower, ymax = upper), position = position_dodge(.50)) + 
  facet_grid(Drug~Sex) +
  scale_col_drug + 
  ylim(0,126) +
  theme(panel.spacing = unit(2, "lines")) + 
  ylab("Dose obtained for the 2nd drug administration")
Obtained dosage by drug, state induction and sex.

Figure 5.1: Obtained dosage by drug, state induction and sex.

my_data_beh %>% 
  .[, .(participant,Sex, state.condition, effect_obtained, Drug)] %>% 
  dcast(participant + Drug + Sex ~state.condition, value.var = "effect_obtained") %>% 
  .[, stress_effect := stress-control] %>% 
  ggplot(aes(x = Sex, y = stress_effect, color = Sex, fill = Sex)) + 
  scale_col_sex + 
  stat_summary(fun.y = mean,alpha = .75,
               geom = "bar", color = NA) + 
  stat_summary(fun.data = mean_cl_normal, color = "black",
               geom = "errorbar", width = .25) + 
  facet_wrap(~Drug) + 
  geom_quasirandom(shape = 1, color = "black") + 
  theme(legend.position = "right")
Effect of stress on obtained effect by sex and drug in raw data.

Figure 5.2: Effect of stress on obtained effect by sex and drug in raw data.

my_data_beh %>% 
  .[, .(participant,Sex, state.condition, effect_obtained, Drug)] %>% 
  dcast(participant + Drug + Sex ~state.condition, value.var = "effect_obtained") %>% 
  .[, stress_effect := stress-control] %>% 
  dcast(participant + Sex ~ Drug, value.var = "stress_effect") %>% 
  .[, stress_effect_x_drug := oxycodone - placebo] %>% 
  ggplot(aes(x = Sex, y = stress_effect_x_drug, color = Sex, fill = Sex)) + 
  scale_col_sex + 
  stat_summary(fun.y = mean,alpha = .75,
               geom = "bar", color = NA) + 
  stat_summary(fun.data = mean_cl_normal, color = "black",
               geom = "errorbar", width = .25) + 
  geom_quasirandom(shape = 1, color = "black") + 
  theme(legend.position = "right")
Interaction of the effct of stress on obtained effect with drug by sex. Larger values indicate greater drug-induced wanting increase in the ocyxodone condition.

Figure 5.3: Interaction of the effct of stress on obtained effect with drug by sex. Larger values indicate greater drug-induced wanting increase in the ocyxodone condition.

Raw data: Select only the oxycodone condition (According to pre-registered plan):

my_data_beh.all = 
  my_data_beh[, .(Sex, participant, session,
                  Drug, state.condition, effect_obtained, target,
                  mean.delta_effect_obtained,
                  first.time.target,n.presses)] 
my_data_beh = 
  copy(my_data_beh.all) %>% 
  .[Drug == "oxycodone"]

To deal with the multi-modal nature of the data, we use an ordered logistic regression, for which we coarsen the data in 25 equally sized ordered categories.

n_cat = 25
my_data_beh = 
  coarsen(my_data_beh,
          var = "effect_obtained",
          bins = n_cat,
          range = c(0,125))

5.1.2 Linear regression without sex effects

We implement a hierarchical regression model where random effects account for repeated measurements. Even though the raw data show that the effect of stress is stronger in men, we first analyse the data without modelling sex effects. The regression model uses fixed effects parameters to model the effects of stress manipulation and random effects to model varying intercepts and stress effects across participants and effects of session number.

my_data_beh[, session := factor(session)]

We estimate a standard hierarchical linear regression with Gaussian errors.

ffit = lme4::lmer(
  effect_obtained ~ state.condition + (1 | participant) + (1|session),
  data = my_data_beh)

To check the model, we first investigate model predictions. The following figure shows that the distribution of model predictions deviates strongly from the observed values and also include impossible negative values.

error_sd = 
  VarCorr(ffit) %>% attr("sc")
p1 = 
  data.table(yhat = predict(ffit)) %>% 
  .[, lower := yhat - error_sd] %>% 
  .[, upper := yhat + error_sd] %>% 
  .[, state.condition := my_data_beh$state.condition] %>% 
  .[order(yhat)] %>% 
  .[, j := 1:.N] %>% 
  ggplot(aes(x = j, y = yhat, ymin = lower, ymax = upper,
             color = state.condition, fill = state.condition)) + 
   geom_hline(yintercept = 0, color = "red") +
  geom_linerange(alpha = .5) +
  geom_point(size = .4) +
  xlab("") +
  ylab("fitted + error") + 
  scale_col_stress
p2 = 
  rbind(my_data_beh[, .(effect_obtained)][, data := "observed"],
           data.table(effect_obtained = predict(ffit), data = "predicted")) %>% 
  ggplot(aes(x = effect_obtained, fill = data)) + 
  geom_histogram(position = "identity", alpha = .5)
  
(p1 | p2) + plot_layout(widths = c(2,1))
Predictions of the linear model +/- 1sd of the normal error distribution.

Figure 5.4: Predictions of the linear model +/- 1sd of the normal error distribution.

Consistently, the following figure shows that the residuals are correlated with the outcome variable and not normally distributed.

res = data.table(residuals = residuals(ffit), fitted = predict(ffit))
p1 = res %>% 
  ggplot(aes(x = fitted, y = residuals)) + 
  geom_point() + 
  geom_smooth()
p2 = 
  res %>% 
  ggplot(aes(x = residuals)) + 
  geom_histogram(aes(y =..density..),
                 breaks = seq(-60, 60, by = 5), 
                 colour = "black", 
                 fill = "white") +
  stat_function(fun = dnorm,
                args = list(mean = mean(res$residuals),
                            sd = sd(res$residuals)))

p1 | p2
Residuals from a linear regression model. Left: Predicted values and residuals are not independent. Right: Residuals are not normally distributed.

Figure 5.5: Residuals from a linear regression model. Left: Predicted values and residuals are not independent. Right: Residuals are not normally distributed.

rm(error_sd,ffit,p1,p2,res)

5.1.3 Ordinal logistic regression (OLR) without sex effects

As for the questionnaire data a linear regression model with Gaussian errors is not a good model for the obtained dosage data. We therefore use ordinal regression models for coarsened data. Specifically, we coarsen the obtained dosages data into 20 equally sized bins and perform Bayesian ordinal logistic regressions on the coarsened data. To report results, we back-calculate effects (including credible intervals) of independent variables on the original 0-125% scale.

We implement a hierarchical ordinal regression model where the cumulative logistic models accounts for the multi-modal distribution of the obtained drug effects and random effects account for repeated measurements. Even though the raw data show that the effect of stress is stronger in men, we first analyse the data without modelling sex effects. The regression model uses fixed effects parameters to model the effects of stress manipulation and random effects to model varying intercepts and stress effects across participants and effects of session number.

analysis_formula = 
  formula(effect_obtained.c ~ state.condition + (1 + state.condition | participant) + (1|session))

We use brms’ default priors for standard deviations of random effects (student-t with 3 df and a sd of 2.5) and a weakly informative normal prior for main effects (normal distribution with mean 0 and standard deviation 2). Priors for the thresholds are specified as normal distributions with the logistic quantile function at cumulative probabilities of response categories as means and standard deviations of 1.

fn = "brmsP1/obtained/nosex.Rdata"
if (file.exists(fn)) {
  load(fn)
} else {
  my_prior = add_cumthresh_prior(
    prior(normal(0,2), class = b),
    my_data_beh$effect_obtained.c)
  fit = brm(
    analysis_formula,
    family = cumulative(),
    data = my_data_beh,
    prior = my_prior,
    backend = "cmdstanr",
    iter = 2000)
  save(fit,file = fn)
}
Click to show sample diagnostics
Table 5.1: Summary of model paramters
parameter mean sd q5 q95 rhat ess_bulk ess_tail
b_Intercept[1] -3.0 0.38 -3.6 -2.4 1.002 4455 2550
b_Intercept[2] -2.6 0.35 -3.2 -2.0 1.000 5015 3363
b_Intercept[3] -2.3 0.32 -2.9 -1.8 1.000 4951 2721
b_Intercept[4] -2.2 0.32 -2.7 -1.7 1.000 5225 2861
b_Intercept[5] -2.1 0.31 -2.6 -1.6 1.000 5232 2949
b_Intercept[6] -1.9 0.30 -2.4 -1.5 1.001 5444 2772
b_Intercept[7] -1.8 0.30 -2.3 -1.3 1.001 5484 2800
b_Intercept[8] -1.5 0.29 -2.0 -1.1 1.000 5556 3330
b_Intercept[9] -1.4 0.29 -1.9 -1.0 1.000 5490 3226
b_Intercept[10] -1.3 0.29 -1.8 -0.9 1.000 5770 3017
b_Intercept[11] -1.2 0.29 -1.6 -0.7 1.000 5282 2919
b_Intercept[12] -1.0 0.29 -1.5 -0.5 1.000 5522 3078
b_Intercept[13] -0.8 0.29 -1.3 -0.3 1.001 5332 3133
b_Intercept[14] -0.6 0.29 -1.1 -0.2 1.000 5555 3412
b_Intercept[15] -0.4 0.30 -0.9 0.1 1.001 5290 3649
b_Intercept[16] -0.1 0.30 -0.6 0.4 1.001 4803 3383
b_Intercept[17] 0.7 0.33 0.2 1.2 1.000 3862 3159
b_Intercept[18] 1.1 0.35 0.5 1.7 1.000 3468 2951
b_Intercept[19] 1.4 0.36 0.8 2.0 1.000 3172 2536
b_Intercept[20] 1.5 0.38 0.9 2.2 1.000 2955 2102
b_state.conditionstress 0.5 0.36 -0.1 1.1 1.001 4779 2654
sd_participant__Intercept 2.1 0.46 1.4 2.9 1.003 930 1523
sd_participant__state.conditionstress 0.7 0.56 0.1 1.8 1.003 793 1007
sd_session__Intercept 0.8 0.53 0.1 1.7 1.002 1410 1672
model: effect_obtained.c ~ state.condition + (1 + state.condition | participant) + (1 | session)
b_ = fix effect parameters
sd_ standard deviation of hierarchical parameters
ess = effective sample size


The next figure shows that the model describes the observed data well.

pp_check(fit, ndraws = 2000, type = "bars")
Posterior predictive check for obtained effects.

Figure 5.6: Posterior predictive check for obtained effects.

5.1.4 OLR without sex effects, adjacent category model with category specific effects

We repeat the same analysis as before, except that we use an adjacent category ordered logistic regression with category specific effects. Differently than the previous analysis, this model does not assume proportional odds. The advantage of this model is that we can fit the data better and have a higher power to detect effects. The disadvantage is that the model estimates for each category transition separate odds ratios. Therefore, it is not meaningful to report one odds ratio as a model parameter.

Here we estimate the model.

analysis_formula = 
  formula(effect_obtained.c ~ cs(state.condition) + (1 + state.condition | participant) + (1|session))
prior_sds = c(.5, 1, 5, 2)
make_prior_b = function(prior_sd) {
  return(
    eval(parse(text = paste("prior(normal(0,",prior_sd,"), class = b)")))
  )
}
for (prior_sd in prior_sds) {
  fn = paste0("brmsP1/obtained/nosex_acat",prior_sd,".Rdata")
  if (file.exists(fn)) {
    load(fn)
  } else {
    my_prior = 
      get_acat_thresh_priors(
        my_data_beh$effect_obtained.c,
        make_prior_b(prior_sd) + 
          prior(normal(0,2), class = "sd"))
    
    fit = brm(
      analysis_formula,
      family = acat(),
      data = my_data_beh,
      prior = my_prior,
      backend = "cmdstanr",
      iter = 2000)
    save(fit,file = fn)
  }
}
Click to show sample diagnostics
Table 5.2: Summary of model paramters
parameter mean sd q5 q95 rhat ess_bulk ess_tail
Threshold[1] 0.7 1.03 -1.0 2.4 1.001 1834 2639
Threshold[2] -2.0 1.11 -3.9 -0.2 1.000 2012 2413
Threshold[3] 0.1 1.20 -1.8 2.2 1.001 3144 2679
Threshold[4] -1.3 1.34 -3.6 0.9 1.003 3554 2433
Threshold[5] -0.9 1.33 -3.1 1.3 1.002 2963 2455
Threshold[6] -2.0 1.25 -4.1 -0.1 1.000 3149 2710
Threshold[7] -1.8 1.05 -3.6 -0.2 1.002 2764 2499
Threshold[8] 0.5 1.10 -1.2 2.4 1.000 3333 2646
Threshold[9] -0.7 1.24 -2.8 1.3 1.002 3640 2867
Threshold[10] -1.7 1.11 -3.6 0.0 1.001 3569 2704
Threshold[11] -0.2 0.96 -1.8 1.4 1.001 3173 2532
Threshold[12] -0.7 0.95 -2.2 0.9 1.000 3067 2509
Threshold[13] -0.5 0.90 -2.0 1.0 1.001 2997 2803
Threshold[14] -0.5 0.83 -1.9 0.9 1.001 3224 2590
Threshold[15] 0.8 0.82 -0.6 2.1 1.001 2912 2238
Threshold[16] -1.1 0.79 -2.4 0.2 1.001 2140 2377
Threshold[17] 1.2 0.71 0.1 2.4 1.002 1737 2377
Threshold[18] 1.3 0.85 0.0 2.8 1.001 2039 2048
Threshold[19] 2.1 1.04 0.5 3.9 1.001 2007 2476
Threshold[20] -1.4 1.02 -3.1 0.3 1.003 1351 2180
bcs_state.conditionstress[1] 2.2 1.22 0.2 4.3 1.002 2657 2932
bcs_state.conditionstress[2] -0.5 1.33 -2.7 1.7 1.000 3301 2684
bcs_state.conditionstress[3] 1.1 1.45 -1.3 3.4 1.000 3443 2927
bcs_state.conditionstress[4] -0.3 1.55 -2.9 2.2 1.001 3714 2947
bcs_state.conditionstress[5] 1.4 1.56 -1.1 4.0 1.000 3775 2714
bcs_state.conditionstress[6] 0.1 1.43 -2.3 2.4 1.001 3522 2681
bcs_state.conditionstress[7] 0.2 1.29 -1.9 2.3 1.001 3582 2999
bcs_state.conditionstress[8] -0.8 1.48 -3.3 1.5 1.001 4293 2917
bcs_state.conditionstress[9] 0.3 1.59 -2.3 2.9 1.000 4269 2628
bcs_state.conditionstress[10] 1.3 1.49 -1.1 3.8 1.001 4254 3015
bcs_state.conditionstress[11] 1.2 1.27 -0.8 3.3 1.002 3915 3206
bcs_state.conditionstress[12] -0.5 1.26 -2.6 1.6 1.000 3674 3088
bcs_state.conditionstress[13] -0.7 1.34 -2.9 1.5 1.000 4276 3131
bcs_state.conditionstress[14] 1.3 1.34 -0.9 3.5 1.001 3599 3095
bcs_state.conditionstress[15] 2.4 1.11 0.6 4.3 1.001 3367 2842
bcs_state.conditionstress[16] -0.2 0.89 -1.7 1.2 1.000 2910 2595
bcs_state.conditionstress[17] -0.5 0.87 -1.9 1.0 1.001 3826 3131
bcs_state.conditionstress[18] -0.5 1.12 -2.3 1.3 1.001 3200 2915
bcs_state.conditionstress[19] 0.5 1.29 -1.5 2.7 1.001 3402 2874
bcs_state.conditionstress[20] 0.3 1.19 -1.6 2.3 1.001 3009 2847
sd_participant__Intercept 1.6 0.55 0.8 2.6 1.007 528 763
sd_participant__state.conditionstress 1.8 0.59 0.9 2.9 1.004 678 1423
sd_session__Intercept 0.6 0.47 0.1 1.6 1.014 635 930
model: effect_obtained.c ~ cs(state.condition) + (1 + state.condition | participant) + (1 | session)
b_ = fix effect parameters
sd_ standard deviation of hierarchical parameters
ess = effective sample size


Next, we calculate the effects on the original 0-125 scale from the fitted model. This approach removes the effect of session by calculating model-predicted outcomes for the two stress conditions in under the assumption that each condition could be in session 1, 2, 3, or 4.

my_contrasts = plot_results_origscale2(fit)
my_contrasts$plot
Estimated mean drug wanting in the 2 state induction conditions (left) and effect of stress, calculate as drug wanting in the stress condition minus drug wanting in the conrol condition (left)

Figure 5.7: Estimated mean drug wanting in the 2 state induction conditions (left) and effect of stress, calculate as drug wanting in the stress condition minus drug wanting in the conrol condition (left)

my_caption = "Model-estimated effect obtained by condition"
my_contrasts$simple.stats %>% 
  kable(caption = my_caption) %>% 
  kable_styling(full_width = FALSE) 
Table 5.3: Model-estimated effect obtained by condition
state.condition m (CrI)
stress 89 (86, 92)
control 84 (80, 87)

In the stress condition participants worked on average for 5 (CI = [1, 10], P_E>0 = 0.99) percentage points higher oxycodone dosage.

5.1.5 Prior sensitivity analysis

To investigate the effect of the model priors on the results we compare the effects calculated with the posterior distribution of the parameters with effects calculate from the prior distribution of the parameters.

First, we generate samples from the prior distribution only and we calculate the same contrasts that we previously calculated with the posterior distributions, this time using the prior distribution.

psa_contrasts = c()
for (prior_sd in prior_sds) {
  fn = paste0("brmsP1/obtained/nosex_acat",prior_sd,"_psa.Rdata")
  if (file.exists(fn)) {
    load(fn)
  } else {
    my_prior = 
      eval(parse(text = paste("prior(normal(0,",prior_sd,"), class = b)")))+ 
      prior(normal(0,1), class = Intercept)
    prior_samples = brm(
      analysis_formula,
      sample_prior = "only", # only prior samples
      family = acat(),
      data = my_data_beh,
      prior = my_prior,
      backend = "cmdstanr",
      iter = 5000)
    save(prior_samples,file = fn)
  }
  # calculate contrasts
  my_contrasts_prior = plot_results_origscale2(prior_samples)
  load(paste0("brmsP1/obtained/nosex_acat",prior_sd,".Rdata"))
  my_contrasts = plot_results_origscale2(fit)
  
  # collect all samples
  psa_contrasts = rbind(
    psa_contrasts,
    rbind(
      my_contrasts$dt.samples.d[, .(value)][, from := "posterior"][, prior_sd := prior_sd],
      my_contrasts_prior$dt.samples.d[, .(value)][, from := "prior"][, prior_sd := prior_sd]
    )
  )
}

And we compare stress effects from posterior and prior under the chosen prior of \(\textrm{Normal}(0,2)\) for the effect of state condition.

psa_contrasts %>% 
  .[, from := factor(from, levels = c("prior","posterior"))] %>% 
  .[prior_sd == 2] %>% 
  ggplot(aes(x = value, color = from)) + 
  geom_density() + 
  xlab("calculated effect") +
  ylab("") +
  theme(legend.position = "right") +
  gg_no_y_axis
Effect distribution under posterior samples and prior prior samples.

Figure 5.8: Effect distribution under posterior samples and prior prior samples.

The figure above illustrates that the prior distribution for the effect is a shrinkage prior that favors null effects over alternative effects. In comparison, the posterior distribution of the effect is narrower and has a non-zero mean. Generally, the employed shrinkage prior means that our estimate of the stress effect is conservative.

Here are calculated effects given the prior and posterior distribution for alternative standard deviations for the prior on the state-condition effect:

psa_contrasts %>% 
  .[, `prior sd` := factor(prior_sd, levels = c("2","0.5","1","5"))] %>% 
  ggplot(aes(x = value, color = from, lty = `prior sd`)) + 
  geom_density() + 
  ylab("calculated effect") +
  xlab("") +
  facet_wrap(~from, scales = "free", nrow = 2) +
  theme(legend.position = "right") +
  gg_no_y_axis
Effect distribution under posterior samples and prior prior samples.

Figure 5.9: Effect distribution under posterior samples and prior prior samples.

psa_contrasts %>% 
  .[, .(`mean (CI)` = get_mci(value,digits = 1)), by = .(from,prior_sd)] %>% 
  dcast(prior_sd ~ from) %>% 
  kable(caption = "Effects of stress condition caculated from prior and posterior distribution.") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic()
Table 5.4: Effects of stress condition caculated from prior and posterior distribution.
prior_sd prior posterior
0.5 -0.3 [-21.3, 19.9]; 0.4849 4.2 [-0.2, 8.7]; 0.97
1.0 -0.3 [-21.8, 20.0]; 0.4896 4.7 [0.1, 9.2]; 0.9777
2.0 0.2 [-23.0, 24.1]; 0.5097 5.1 [0.8, 9.5]; 0.991
5.0 0.8 [-34.5, 35.4]; 0.521 5.7 [1.5, 10.2]; 0.99475

Priors with smaller standard deviations for the effect of stress lead to posterior distributions for the effect of stress with somewhat smaller means and standard deviations, whereas priors with larger standard deviations for the effect of stress would lead to posterior distributions for the effect of stress with somewhat larger mean and standard deviations. In sum, our results show that we have employed shrinkage priors that result in a conservative estimate of the effect of stress.

5.1.6 OLR with sex effects

Because the initial data visualisation also shows sex effects, we estimate sex specific main effects and effects of stress in the next model.

This analysis uses an adjacent category model with category-specific effects, which relaxes the assumption of proportional odds. The basic analysis model is to estimate the the obtained effect dependent on state-condition and Sex, while accounting for repeated measurement within individuals and for session effects with random effects.

In particular we use the following model:

analysis_formula = 
  formula(effect_obtained.c ~ 
            # category specific effects of state condition and sex
            cs(state.condition*Sex) +
            # random effects for participants
            (1 + state.condition | participant) + 
            # random effects for session
            (1|session))
for (prior_sd in prior_sds) {
  fn = paste0("brmsP1/obtained/acat",prior_sd,".Rdata")
  if (file.exists(fn)) {
    load(fn)
  } else {
    my_prior = 
      get_acat_thresh_priors(
        my_data_beh$effect_obtained.c,
        make_prior_b(prior_sd) + 
          prior(normal(0,2), class = "sd"))
    control = list(adapt_delta = .90)
    fit = brm(
      analysis_formula,
      family = acat(),
      data = my_data_beh,
      prior = my_prior,
      backend = "cmdstanr",
      iter = 2000)
    save(fit,file = fn)
    prior_samples = update(fit,sample_prior = "only")
    save(prior_samples,file = gsub(".Rdata","_psa.Rdata",fn))
  }
}
load(paste0("brmsP1/obtained/acat",2,".Rdata"))
Click to show sample diagnostics
Table 5.5: Summary of model paramters
parameter mean sd q5 q95 rhat ess_bulk ess_tail
Threshold[1] 0.9 1.17 -1.0 2.9 1.001 2812 1945
Threshold[2] -2.4 1.26 -4.5 -0.4 1.003 2357 2792
Threshold[3] 0.7 1.35 -1.5 2.9 1.000 2954 2762
Threshold[4] -1.0 1.45 -3.4 1.3 1.000 2901 3177
Threshold[5] -0.9 1.49 -3.4 1.6 1.001 3218 2476
Threshold[6] -2.2 1.40 -4.5 0.1 1.000 3216 2782
Threshold[7] -2.4 1.28 -4.6 -0.4 1.000 3112 3274
Threshold[8] 0.9 1.29 -1.2 3.0 1.001 3164 2542
Threshold[9] -0.6 1.36 -2.8 1.7 1.001 2556 2173
Threshold[10] -2.0 1.29 -4.1 0.1 1.002 2900 2677
Threshold[11] -0.2 1.19 -2.1 1.8 1.001 2680 2677
Threshold[12] -1.2 1.19 -3.2 0.7 1.001 2702 2501
Threshold[13] 0.1 1.11 -1.7 1.9 1.002 2524 2963
Threshold[14] -1.0 1.08 -2.8 0.8 1.003 2435 2681
Threshold[15] 0.8 1.04 -0.9 2.5 1.003 2244 2611
Threshold[16] -1.4 0.97 -3.1 0.1 1.001 2010 2748
Threshold[17] 0.9 0.85 -0.5 2.3 1.004 1820 2086
Threshold[18] 1.2 1.01 -0.4 2.9 1.001 1875 2556
Threshold[19] 2.2 1.17 0.4 4.2 1.001 2026 2038
Threshold[20] -2.2 1.15 -4.1 -0.3 1.003 1481 2029
bcs_state.conditionstress[1] 1.3 1.38 -1.0 3.6 1.001 2661 2492
bcs_state.conditionstress[2] -0.1 1.47 -2.5 2.4 1.000 3571 2944
bcs_state.conditionstress[3] 0.2 1.59 -2.5 2.8 1.000 3419 2726
bcs_state.conditionstress[4] -0.8 1.67 -3.6 1.9 1.001 3613 2578
bcs_state.conditionstress[5] 0.8 1.60 -1.9 3.4 1.001 3534 2640
bcs_state.conditionstress[6] -0.2 1.63 -3.0 2.5 1.001 3833 2801
bcs_state.conditionstress[7] 0.2 1.56 -2.4 2.7 1.002 3205 2992
bcs_state.conditionstress[8] -0.7 1.65 -3.5 2.0 1.000 4074 3000
bcs_state.conditionstress[9] 0.2 1.67 -2.5 3.0 1.004 4024 2986
bcs_state.conditionstress[10] 1.1 1.69 -1.6 3.9 1.000 4231 2695
bcs_state.conditionstress[11] 1.1 1.57 -1.4 3.7 1.000 3540 2658
bcs_state.conditionstress[12] 0.3 1.51 -2.1 2.8 1.000 3193 2998
bcs_state.conditionstress[13] 0.0 1.55 -2.5 2.5 1.001 3452 2586
bcs_state.conditionstress[14] 1.8 1.54 -0.8 4.3 1.002 3642 3019
bcs_state.conditionstress[15] 2.9 1.34 0.7 5.1 1.001 3345 2610
bcs_state.conditionstress[16] 0.9 1.10 -0.9 2.7 1.002 2714 2697
bcs_state.conditionstress[17] -0.2 1.01 -1.9 1.4 1.001 2480 2407
bcs_state.conditionstress[18] -0.8 1.24 -2.8 1.3 1.001 2381 2393
bcs_state.conditionstress[19] 0.5 1.36 -1.7 2.7 1.000 3201 2689
bcs_state.conditionstress[20] 0.8 1.31 -1.4 2.9 1.002 2719 2803
bcs_Sexwoman[1] 0.5 1.28 -1.7 2.5 1.001 2753 2532
bcs_Sexwoman[2] -0.2 1.39 -2.5 2.1 1.000 2940 3041
bcs_Sexwoman[3] 1.2 1.51 -1.3 3.7 1.002 3030 2175
bcs_Sexwoman[4] 0.9 1.53 -1.6 3.4 1.002 2870 2645
bcs_Sexwoman[5] 0.0 1.58 -2.6 2.5 1.002 3070 2666
bcs_Sexwoman[6] 0.1 1.49 -2.4 2.5 1.000 3428 3054
bcs_Sexwoman[7] -1.0 1.41 -3.3 1.3 1.001 2910 2964
bcs_Sexwoman[8] 1.1 1.51 -1.4 3.6 1.001 3183 2893
bcs_Sexwoman[9] 0.3 1.51 -2.2 2.8 1.001 2925 2835
bcs_Sexwoman[10] -0.5 1.46 -2.9 1.9 1.002 2482 2455
bcs_Sexwoman[11] 0.0 1.34 -2.2 2.2 1.002 2953 2868
bcs_Sexwoman[12] -1.2 1.36 -3.4 1.0 1.001 2923 2849
bcs_Sexwoman[13] 1.2 1.29 -0.8 3.4 1.001 2524 2911
bcs_Sexwoman[14] -1.1 1.22 -3.1 0.9 1.003 2678 2612
bcs_Sexwoman[15] -0.3 1.17 -2.1 1.7 1.001 2503 3028
bcs_Sexwoman[16] -0.8 1.11 -2.6 1.0 1.002 2425 2668
bcs_Sexwoman[17] -1.1 1.03 -2.9 0.6 1.002 2200 2769
bcs_Sexwoman[18] -1.0 1.18 -3.0 0.9 1.001 2638 2895
bcs_Sexwoman[19] -0.7 1.29 -2.8 1.4 1.001 2069 2778
bcs_Sexwoman[20] -2.7 1.25 -4.8 -0.7 1.001 2394 2678
bcs_state.conditionstress:Sexwoman[1] 1.2 1.44 -1.1 3.6 1.001 2818 2857
bcs_state.conditionstress:Sexwoman[2] -1.5 1.61 -4.2 1.0 1.000 3049 2846
bcs_state.conditionstress:Sexwoman[3] 1.2 1.66 -1.6 3.9 1.000 4037 3012
bcs_state.conditionstress:Sexwoman[4] -0.2 1.66 -3.0 2.4 1.000 3621 2821
bcs_state.conditionstress:Sexwoman[5] 1.1 1.71 -1.7 3.9 1.001 3247 2751
bcs_state.conditionstress:Sexwoman[6] -0.1 1.64 -2.7 2.6 1.002 3353 2489
bcs_state.conditionstress:Sexwoman[7] 0.2 1.60 -2.4 2.8 1.000 3382 2712
bcs_state.conditionstress:Sexwoman[8] -1.3 1.69 -4.1 1.4 1.001 3683 2905
bcs_state.conditionstress:Sexwoman[9] -0.4 1.72 -3.3 2.4 1.002 3972 2711
bcs_state.conditionstress:Sexwoman[10] 0.5 1.72 -2.4 3.3 1.001 4229 3165
bcs_state.conditionstress:Sexwoman[11] 0.2 1.60 -2.5 2.8 1.003 3218 2645
bcs_state.conditionstress:Sexwoman[12] -1.0 1.53 -3.5 1.5 1.000 3321 2972
bcs_state.conditionstress:Sexwoman[13] -1.9 1.60 -4.5 0.7 0.999 3896 2985
bcs_state.conditionstress:Sexwoman[14] -0.3 1.62 -3.0 2.4 1.000 4040 3062
bcs_state.conditionstress:Sexwoman[15] -0.4 1.48 -2.9 2.0 1.002 3331 3084
bcs_state.conditionstress:Sexwoman[16] -2.1 1.28 -4.2 0.0 1.000 2797 2890
bcs_state.conditionstress:Sexwoman[17] -1.3 1.34 -3.5 0.8 1.002 3194 2663
bcs_state.conditionstress:Sexwoman[18] 1.3 1.56 -1.3 3.9 1.000 3340 2396
bcs_state.conditionstress:Sexwoman[19] 0.0 1.58 -2.5 2.7 1.001 3533 2876
bcs_state.conditionstress:Sexwoman[20] -1.0 1.46 -3.4 1.4 1.001 3604 3202
sd_participant__Intercept 1.9 0.59 1.1 3.0 1.005 406 1003
sd_participant__state.conditionstress 2.0 0.67 1.1 3.2 1.009 413 719
sd_session__Intercept 0.9 0.60 0.2 2.0 1.003 486 725
model: effect_obtained.c ~ cs(state.condition * Sex) + (1 + state.condition | participant) + (1 | session)
b_ = fix effect parameters
sd_ standard deviation of hierarchical parameters
ess = effective sample size


To communicate the results we calculate conditional effects, which show by how many percentage points the obtained effects increased on average in the stress condition. See the function make_plot_samples in the utils.R file for details.

plot.stats_table_obtained = 
  plot_results_origscale(
    fit,
    my_data_beh,
    outcome.var = "effect_obtained.c")

p1 = plot.stats_table_obtained$p1
yrange = layer_scales(p1)$y$range$range
p1 = p1 + scale_y_continuous(limits = yrange) +
  guides(y = guide_axis_truncated(
    trunc_lower = yrange[1]+diff(yrange)/500)) +
  coord_cartesian(ylim = c(67,108.5)) +
  xlab("") +
  theme(plot.title = element_text(hjust = -.05))

p2 = plot.stats_table_obtained$p2

px = 
  (p1 | p2) + 
  plot_layout(widths = c(4,1),) +
  plot_annotation(tag_levels = 'A')
ggsave(px,
        filename = "figures/Effect_Obtained.png", 
        width = 20,height = 15, units = "cm")
ggsave(px,
        filename = "figures/Effect_Obtained.pdf", 
        width = 20,height = 15, units = "cm", device = cairo_pdf)
“Effect of stress on effect obtained in the self administration task”
“Effect of stress on effect obtained in the self administration task”
my_caption = "Model-estimated effect obtained by sex and state condition."
my_footnote = "Numbers are means and CIs: Upper and lower bound of 95% credible intervals"
plot.stats_table_obtained$simple.stats %>% 
  kable(digits = 2, caption = my_caption) %>% 
  kable_styling(full_width = FALSE) %>%
  kable_classic() %>% 
  add_footnote(my_footnote, notation="none")
Table 5.6: Model-estimated effect obtained by sex and state condition.
Sex control stress
men 91 (86, 96) 104 (99, 108)
total 84 (81, 87) 89 (86, 92)
women 77 (72, 81) 74 (69, 78)
Numbers are means and CIs: Upper and lower bound of 95% credible intervals
plot.stats_table_obtained$tbl_me_sex
Table 5.7: Effect of sex, calculated as effect obtained in men - women.
state condition effect_woman
control 14 [8, 20]; >.999
stress 30 [24, 37]; >.999
overall 22 [18, 27]; >.999
Numbers are means and CIs: Upper and lower bound of 95% credible intervals
my_caption = "Effect of stress on effect obtained by sex, over all, and sex difference"
my_footnote = "CI: Upper and lower bound of 95% credible intervals. P > 0: Posterior probability that the parameter or contrast is larger than 0."
plot.stats_table_obtained$stats[,.(Sex, mean, CI, `P > 0`)] %>% 
  .[, Sex := factor(Sex, levels = c("total","women","men","sex diff: m-f"))] %>% 
  setkeyv("Sex") %>% 
  kable(digits = 2, caption = my_caption) %>% 
  kable_styling(full_width = FALSE) %>%
  kable_classic() %>% 
  add_footnote(my_footnote, notation="none")
Table 5.8: Effect of stress on effect obtained by sex, over all, and sex difference
Sex mean CI P > 0
total 4.91 0.56, 9.22 0.99
women -3.07 -9.1, 3.33 0.16
men 12.89 6.72, 18.66 1.00
sex diff: m-f 15.96 7.02, 24.49 1.00
CI: Upper and lower bound of 95% credible intervals. P > 0: Posterior probability that the parameter or contrast is larger than 0.
PS = plot.stats_table_obtained$stats

We used hierarchical ordinal logistic regression models to estimate the effect of stress on the drug wanting, which was measured as the percent of the drug effect from the first drug administration participants worked for in the self administration task. To better understand the magnitude of these effects, we calculated marginal effects of stress on drug wanting. Taking sex differences into account revealed that whereas stress over all increased wanting by 4.9 (CI=[0.56, 9.22], BF=69.2) points, stress increased drug wanting in only males. In males, wanting was 12.9 (CI=[6.72, 18.66], BF=Inf) percentage points higher in the stress compared to the control condition, and in females the differences was -3.1 (CI=[-9.1, 3.33], BF=0.2). The male-female differences of 16 (CI=[7.02, 24.49], BF=Inf) percentage points indicates that the effect of stress on drug wanting was higher in males than in females.

To investigate potential explanations for the sex difference in drug wanting we explored sex differences in the effect of drug on positive and negative affect, embarrassment/selfconsciousness and drug-side. This exploratory analysis showed that while females reported greater side effects (no differences for other affect variables), including drug effects on these potential mediators did not abolish the sex differences. This suggests that sex differences of drug effects on affect and side effects cannot fully explain the sex differences in drug wanting.

5.2 Effect of drug on post-stress induction affect

The following plots show changes in affect from pre drug administration (stage 4) to post drug administration (stage 5). The effect of drug is encoded in different pre-post slopes for the placebo and oxycodone conditions and can very by sex.

my_data.QPS = 
  my_data.Q %>% 
  .[, Stage := factor(stage.names[stage], levels = stage.names)] %>% 
  .[stage %in% c(4,5) & 
              item_name %in% c(negaffect_items,
                               posaffect_items,
                               embselfcn_items)] %>% 
  .[, affect := ifelse(item_name %in% negaffect_items,
                       "negative",
                       ifelse(item_name %in% posaffect_items,
                              "positive","emb_selfc"))] %>% 
  coarsen(var = "response",
          bins = 20,
          range = c(0,100))
my_data.QPS_stats = 
  my_data.QPS %>% 
  .[, .(response = mean(response),
        se = sd(response)/sqrt(.N)),
    by = .(Drug, stage, affect, Sex)] %>% 
  .[, `:=`(lower = response-2*se,
           upper = response+2*se)] 

my_data.QPS_stats %>% 
  ggplot(aes(x = stage, y = response, color = Drug,
             ymin = lower, ymax = upper, shape = Sex, lty = Sex)) +
  geom_point(position = pd.2) + 
  geom_line(aes(group = factor(Drug):factor(Sex)), position = pd.2) +
  geom_linerange(position = pd.2) + 
  facet_wrap(~affect) + 
  scale_x_continuous(breaks = c(3.75,4.75), labels = c("state\ninduction","dose1"), limits = c(3.5,5.25)) +
  scale_col_drug

The means already show clearly affect did not develop more positively in the oxycodone condition compared to the placebo condition.

We can make the same comparison for side effects:

my_data.QSE = 
  my_data.Q %>% 
  .[stage %in% c(4,5) & 
              item_name %in% sideffect_items] %>% 
  coarsen(var = "response",
          bins = 20,
          range = c(0,100))

my_data.QSE_stats = 
  my_data.QSE %>% 
  .[, .(response = mean(response),
        se = sd(response)/sqrt(.N)),
    by = .(Drug, stage, item_name, Sex)] %>% 
  .[, `:=`(lower = response-2*se,
           upper = response+2*se)] 

my_data.QSE_stats %>% 
    ggplot(aes(x = stage, y = response, color = Drug,
               ymin = lower, ymax = upper, shape = Sex, lty = Sex)) +
    geom_point(position = pd.2) + 
    geom_line(aes(group = factor(Drug):factor(Sex)), position = pd.2) +
    geom_linerange(position = pd.2) + 
    facet_wrap(~item_name) + 
  scale_x_continuous(breaks = c(3.75,4.75), labels = c("state\ninduction","dose1"), limits = c(3.5,5.25)) +
  scale_col_drug

This figure suggests that there was a greater increase in negative drug-related emotions in women than in men.

5.3 OLR with sex effects for placebo condition

my_data_beh_placebo = 
  copy(my_data_beh.all) %>% 
  .[Drug == "placebo"] %>% 
  .[effect_obtained > 125, effect_obtained := 125] %>% 
  coarsen(
    var = "effect_obtained",
    bins = n_cat,
    range = c(0,125))
analysis_formula = 
  formula(effect_obtained.c ~ state.condition:Sex + Sex + (1 + state.condition | participant) + (1|session))

We fit the model for effects obtained

options(mc.cores = 4)
fn = "brmsP1/obtained/placebo.Rdata"
if (file.exists(fn)) {
  load(fn)
} else {
  my_prior = 
    get_acat_thresh_priors(
      my_data_beh_placebo$effect_obtained.c,
      make_prior_b(prior_sd) + 
        prior(normal(0,2), class = "sd"))
  fit = brm(
    analysis_formula,
    family = acat(),
    data = my_data_beh_placebo,
    prior = my_prior,
    backend = "cmdstanr",
    iter = 2000)
  save(fit,file = fn)
}

And we make a figure as we did for the oxycodone condition:

plot.stats_table_obtained_placebo = 
  plot_results_origscale(
    fit,
    my_data_beh_placebo,
    outcome.var = "effect_obtained.c")

ggsave(plot.stats_table_obtained_placebo$plot,
       filename = "figures/EOplac.png", 
       width = 15,height = 15, units = "cm", device = cairo_pdf)

5.4 Stated desired effect

Prepare data

my_data_beh = 
  merge(my_data_beh,
        cheaters,
        by = c("participant","Drug","state.condition","session"),
        all.x = TRUE)
my_data_beh = 
  coarsen(my_data_beh,
          var = "target",
          bins = n_cat,
          range = c(0,125))

my_data_beh[is.na(cheater)] %>% 
  ggplot(aes(x = target, fill = state.condition)) + 
  geom_histogram(alpha = .5, position = "identity") + 
  scale_col_stress + 
  facet_grid(Sex~Drug) 

analysis_formula = 
  formula(target.c ~ cs(state.condition*Sex) + (1 + state.condition | participant) + (1|session))
fn = "brmsP1/desired.Rdata"
if (file.exists(fn)) {
  load(fn)
} else {
  my_prior = 
      prior(normal(0,2), class = b) + 
      prior(normal(0,1), class = Intercept)
  fit = brm(
    analysis_formula,
    family = acat(),
    data = my_data_beh[is.na(cheater)],
    prior = my_prior,
    iter = 5000,
    backend = "cmdstanr")
  save(fit,file = fn)
}

Make the plots

plot.stats_table_target = 
  plot_results_origscale(
    fit,
    my_data_beh[is.na(cheater)],
    outcome.var = "target.c")

p1 = plot.stats_table_target$p1 + ylab("% desired effect (0-125)")
yrange = layer_scales(p1)$y$range$range
p1 = p1 + scale_y_continuous(limits = yrange) +
  guides(y = guide_axis_truncated(
    trunc_lower = yrange[1]+diff(yrange)/500)) +
  coord_cartesian(ylim = c(67,108.5)) +
  xlab("") +
  theme(plot.title = element_text(hjust = -.05))

p2 = plot.stats_table_target$p2

px = 
  (p1 | p2) + 
  plot_layout(widths = c(4,1),) +
  plot_annotation(tag_levels = 'A')
ggsave(px,
        filename = "figures/Effect_Desired.png", 
        width = 20,height = 15, units = "cm")
ggsave(px,
        filename = "figures/Effect_Desired.pdf", 
        width = 20,height = 15, units = "cm")

Here are the estimated responses:

my_caption = "Model-estimated stated desired effect by sex and state condition."
my_footnote = "Numbers are means and CIs: Upper and lower bound of 95% credible intervals"
plot.stats_table_target$simple.stats %>% 
  kable(digits = 2, caption = my_caption) %>% 
  kable_styling() %>% 
  add_footnote(my_footnote, notation="none")
Table 5.9: Model-estimated stated desired effect by sex and state condition.
Sex control stress
men 90 (85, 95) 100 (95, 104)
total 84 (80, 88) 88 (85, 91)
women 78 (74, 83) 77 (72, 81)
Numbers are means and CIs: Upper and lower bound of 95% credible intervals

Here are the estimated effects:

“Effect of stress on desired effect in the self administration task”
“Effect of stress on desired effect in the self administration task”

5.5 Take drug again

Prepare data

my_data.Q.take_again = 
  my_data.Q %>% 
  .[Drug == "oxycodone" & item_name == "drug_again"] %>% 
  .[, .(participant,session,Sex,Drug,state.condition,stage, response)] %>% 
  coarsen(var = "response",
          bins = 20,
          range = c(0,100))

my_data.Q %>% 
  .[item_name == "drug_again"] %>% 
  ggplot(aes(x = response, fill = state.condition)) + 
  geom_histogram(alpha = .5, position = "identity") + 
  scale_col_stress + 
  facet_grid(Sex~Drug)
Willingess to take the drug again by gender and state condition.

Figure 5.10: Willingess to take the drug again by gender and state condition.

Estimate the model (only with data from the oxycodon condition).

analysis_formula = 
  formula(response.c ~ cs(state.condition*Sex) + (1 + state.condition | participant) + (1|session))
options(mc.cores = 4)

fn = "brmsP1/again.Rdata"
if (file.exists(fn)) {
  load(fn)
} else {
  fit = brm(
    analysis_formula,
    family = acat(),
    data = my_data.Q.take_again,
    prior = my_prior,
    iter = 5000,
    backend = "cmdstanr")
  save(fit,file = fn)
}

Make the plots

plot.stats_table_again = 
  plot_results_origscale(
    fit,
    my_data.Q.take_again,
    outcome.var = "response.c")

Here are the estimated responses:

my_caption = "Model-estimated 'take again' responses by sex and state condition."
my_footnote = "Numbers are means and CIs: Upper and lower bound of 95% credible intervals"
plot.stats_table_again$simple.stats %>% 
  kable(digits = 2, caption = my_caption) %>% 
  kable_styling() %>% 
  add_footnote(my_footnote, notation="none")
Table 5.10: Model-estimated ‘take again’ responses by sex and state condition.
Sex control stress
men 31 (29, 33) 30 (27, 33)
total 28 (27, 30) 29 (27, 32)
women 26 (23, 29) 28 (26, 30)
Numbers are means and CIs: Upper and lower bound of 95% credible intervals

Here are the estimated effects:

cpt = "Effects of stress on 'take drug again'"
plot.stats_table_again$stats %>% 
  .[, CI := paste0("(",round(lower),", ", round(upper),")")] %>% 
  .[, c(1,3,8,7,6)] %>% 
  .[, `P > 0` := round(`P > 0`,2)] %>%
  .[, mean := round(mean)] %>%
  kable(caption = cpt) %>% 
  kable_styling(full_width = FALSE)
Table 5.11: Effects of stress on ‘take drug again’
Sex mean CI BF P > 0
men -1 (-5, 2) 0.6 0.38
women 2 (-2, 5) 5.7 0.85
total 1 (-2, 3) 2.0 0.67
sex diff: m-f -3 (-8, 2) 0.2 0.18

5.6 Plot for stated dersired effect and take drug again

Put plots together:

p = 
  (plot.stats_table_obtained$p1 + xlab("") + ylim(0,120) +
     scale_alpha(guide = 'none') +
     guides(y = guide_axis_truncated(trunc_lower = 0, trunc_upper = 120)) +
     theme(legend.position = "none", text = element_text(size = 10)) + 
     ylab("Drug self administration obtained (% of 1st)")) |
  (plot.stats_table_target$p1 + xlab("") + ylim(0,120) +
     guides(y = guide_axis_truncated(trunc_lower = 0, trunc_upper = 120)) +
     scale_alpha(guide = 'none') +
     theme(legend.position = "bottom", text = element_text(size = 10)) + 
     ylab("Drug self administration target (% of 1st)")) | 
  (plot.stats_table_again$p1 + ylim(0,100) +
     guides(y = guide_axis_truncated(trunc_lower = 0, trunc_upper = 100)) +
     scale_alpha(guide = 'none') +
     theme(legend.position = "none", text = element_text(size = 10)) + 
     ylab("Drug take again (%)")) +
  plot_layout(guides = "collect") +
  plot_annotation(tag_levels = 'A')
ggsave(p,filename = "figures/MR.png",width = 20,height = 10, units = "cm")
Figure X: Effect of stress on desired effect, obtained effect and rating to “take the drug again”.
Figure X: Effect of stress on desired effect, obtained effect and rating to “take the drug again”.

6 Distribution of the number of button presses:

my_data_beh.all[, session := factor(session)]

my_data_beh.all = 
  merge(my_data_beh.all,
        cheaters,
        by = c("participant","Drug","state.condition","session"),
        all.x = TRUE) %>% 
  .[, n.presses.s := n.presses-mean(n.presses), by = .(participant)]


my_data_beh.all[, `target effect`  := ifelse(target > 50, "target > 50%","target < 50%")]

hist.presses = 
  my_data_beh.all[is.na(cheater)] %>% 
  ggplot(aes(x = n.presses, fill = Drug)) + 
  geom_histogram(position = "identity", alpha = .5, bins = 15) + 
  scale_fill_manual(values = palette()[c(2,4)]) +
  facet_grid(Sex~`target effect`) + 
  xlab("Number of presses") + 
  ylab("Number of participants") + 
  scale_col_drug

hist.time2target = 
  my_data_beh.all[is.na(cheater)] %>% 
  ggplot(aes(x = first.time.target, fill = Drug)) + 
  geom_histogram(position = "identity", alpha = .5, bins = 15) + 
  scale_fill_manual(values = palette()[c(2,4)]) +
   facet_grid(Sex~`target effect`) + 
  xlab("Time to target") + 
  ylab("Number of participants") + 
  scale_col_drug

#hist.presses / hist.time2target

my_data_beh.all %>% 
  .[, .(m = round(mean(n.presses)), sd = round(sd(n.presses))),
   by = .(Drug,state.condition,Sex)] %>% 
  kable(caption = "Number of button presses by druig and stress condition and sex") %>% 
  kable_styling(full_width = FALSE)
Table 6.1: Number of button presses by druig and stress condition and sex
Drug state.condition Sex m sd
oxycodone control women 613 101
oxycodone stress women 622 64
placebo control women 579 178
placebo stress women 615 138
oxycodone control men 639 161
oxycodone stress men 669 163
placebo control men 656 160
placebo stress men 670 135
# my_data_beh.all %>% 
#   ggplot(aes(x = state.condition, y = n.presses, color = Drug)) + 
#   stat_pointinterval(position = position_dodge(.2),  .width = c(0.5, 0.9)) + 
#   scale_color_manual(values = palette()[c(2,4)]) +
#   facet_wrap(~Sex) + 
#   scale_col_drug

7 Stress relief

To test if oxycodone buffers the effect of stress, we compare if the stress ratings increase more from pre to post state induction in the placebo condition compared to the oxycodone condition.

Here we plot the available data points:

my_data.Q.stressed %>% plot_stages()

The last pre state induction and pre drug time point is “BL 3” and the first post state induction and post drug time point is “dose 1”. Hence, we compare if the stress ratings increase more from “BL 3” to “dose 1” in the placebo condition compared to the oxycodone condition.

my_data.Q.stressdrug = select_data("stressed")
my_stages = c(3,5)
the_data = make_the_data(my_data.Q.stressdrug, "stressdrug")[, outcome := "stress"]

p1 = the_data %>% 
  .[, .(response = mean(response)), 
    by = .(participant, session, Drug, state.condition, stage)] %>% 
  ggplot(aes(x = stage, y = response, color = Drug,
             group = factor(participant):factor(session))) + 
  geom_point(alpha = .5) + 
  facet_wrap(~state.condition) + 
  geom_line(alpha = .5) + 
  scale_col_drug 

pdata = 
  the_data %>% 
  # average the pre stages BL1, BL2, BL3
  .[, .(response = mean(response)), 
    by = .(participant, Drug, state.condition, stage, Sex)] %>% 
  # calculate state induction effect as post-pre difference
  dcast(participant  + Drug + state.condition + Sex ~ stage, value.var = "response") %>%
  .[, induction_effect := post.stressdrug - pre.stressdrug] %>% 
  # calculate stress effect as difference of post-pre difference in stress & control conditions
  dcast(participant + Drug + Sex ~ state.condition, value.var = "induction_effect") %>% 
  .[, stress_effect := stress - control]

pdata.stats = 
  seWithin(pdata, "stress_effect",
           withinvars = c("Drug"),
           betweenvars = c("Sex"),
           idvar = "participant") %>% 
  data.table() %>% 
  .[, `:=`(lower=stress_effect-ci,
           upper=stress_effect+ci)]

p2 = pdata.stats %>% 
  ggplot(aes(x = Drug, y = stress_effect, fill = Drug)) + 
  geom_bar(stat = "identity", alpha = .5) + 
  geom_errorbar(aes(ymin = lower, ymax = upper), width = .25, size = 1) + 
  facet_wrap(~Sex) + 
  geom_beeswarm(data = pdata, aes(color = Drug)) +
  scale_col_drug +
  theme(legend.position = "none")

p3 = pdata  %>% 
  # calculate drug x stress interaction effect as difference of stress effects between drug conditions
  dcast(participant + Sex ~ Drug, value.var = "stress_effect") %>% 
  .[, stress_x_drug := oxycodone - placebo] %>% 
  ggplot(aes(x = stress_x_drug, fill=Sex)) + 
  geom_histogram(alpha = .5, position = "identity", bins = 20) +
  scale_col_sex +
  theme(legend.position = "right") 

p1 / (p2 | p3 ) + plot_annotation(tag_levels = 'A')
A: Stress ratings from pre to post state induction & drug administration. B: Effect of stress induction in the raw data measured after drug adminstration. Stress effects are calculated as the change in stress ratings from pre to post stress induction and drug administration. C: Distribution of the interaction of stress effect (pre-post induction) with drug condition (directly calculated from the data).

Figure 7.1: A: Stress ratings from pre to post state induction & drug administration. B: Effect of stress induction in the raw data measured after drug adminstration. Stress effects are calculated as the change in stress ratings from pre to post stress induction and drug administration. C: Distribution of the interaction of stress effect (pre-post induction) with drug condition (directly calculated from the data).

We coarsen the outcome variable in preparation of the ordinal regression model

the_data = the_data %>% 
  coarsen(var = "response", bins = 20, range = c(0,100))

stress_drug_formula = stress_drug_formula = response.c ~ induced.state:Sex:Drug + (1|session) + (1|participant/stage) + (induced.state:Sex:Drug|participant)

stress_drug_formula = stress_drug_formula = response.c ~ induced.state:Sex:Drug + (1|session) + (1|participant/stage) + (0+induced.state:Sex:Drug|participant)

We fit the model

fn = "brmsP1/stress_relief_acat.Rdata"
if (file.exists(fn)) {
  load(fn)
} else {
  fit = fit_ordinal_model(the_data, stress_drug_formula, fn, family = acat())
}
#pp_check(fit, ndraws = 2000, type = "bars")
Click to show sample diagnostics
Table 7.1: Summary of model paramters
parameter mean sd q5 q95 rhat ess_bulk ess_tail
Threshold[1] 1.0 0.73 -0.2 2.2 1.002 784 1312
Threshold[2] 0.5 0.73 -0.7 1.7 1.002 803 1438
Threshold[3] 0.5 0.74 -0.8 1.7 1.002 821 1259
Threshold[4] 0.9 0.75 -0.3 2.2 1.002 839 1484
Threshold[5] 0.5 0.76 -0.8 1.7 1.002 869 1668
Threshold[6] 0.6 0.78 -0.7 1.9 1.002 872 1516
Threshold[7] 0.6 0.79 -0.6 1.9 1.003 934 1559
Threshold[8] 0.9 0.82 -0.4 2.2 1.002 910 1587
Threshold[9] 0.2 0.82 -1.1 1.5 1.002 976 1753
Threshold[10] 1.4 0.86 -0.1 2.8 1.002 1017 1901
Threshold[11] 0.9 0.92 -0.6 2.4 1.002 1146 2017
Threshold[12] 1.2 1.03 -0.5 2.9 1.002 1273 2217
Threshold[13] 0.1 1.00 -1.6 1.7 1.001 1368 2118
Threshold[14] 2.5 1.20 0.6 4.6 1.001 1462 2341
Threshold[15] 0.6 1.39 -1.6 3.0 1.001 2255 2693
b_induced.statecontrol:Sexman:Drugoxycodone -0.1 0.70 -1.3 1.0 1.002 811 1459
b_induced.statestress:Sexman:Drugoxycodone 0.0 0.72 -1.2 1.1 1.002 808 1361
b_induced.statecontrol:Sexwoman:Drugoxycodone 0.0 0.70 -1.2 1.1 1.003 802 1433
b_induced.statestress:Sexwoman:Drugoxycodone 0.2 0.71 -1.0 1.3 1.003 809 1428
b_induced.statecontrol:Sexman:Drugplacebo -0.1 0.71 -1.3 1.0 1.003 799 1454
b_induced.statestress:Sexman:Drugplacebo 0.0 0.71 -1.1 1.2 1.003 814 1438
b_induced.statecontrol:Sexwoman:Drugplacebo -0.1 0.71 -1.2 1.1 1.003 804 1415
b_induced.statestress:Sexwoman:Drugplacebo 0.2 0.71 -0.9 1.4 1.004 791 1459
sd_participant__Intercept 0.3 0.07 0.2 0.4 1.006 826 799
sd_participant__induced.statecontrol:Sexman:Drugoxycodone 0.3 0.15 0.1 0.6 1.001 711 722
sd_participant__induced.statestress:Sexman:Drugoxycodone 0.2 0.18 0.0 0.6 1.003 1489 1962
sd_participant__induced.statecontrol:Sexwoman:Drugoxycodone 0.3 0.12 0.1 0.5 1.002 957 781
sd_participant__induced.statestress:Sexwoman:Drugoxycodone 0.3 0.17 0.0 0.6 1.006 784 1154
sd_participant__induced.statecontrol:Sexman:Drugplacebo 0.2 0.14 0.0 0.5 1.007 578 1194
sd_participant__induced.statestress:Sexman:Drugplacebo 0.2 0.14 0.0 0.5 1.007 1629 2090
sd_participant__induced.statecontrol:Sexwoman:Drugplacebo 0.3 0.14 0.1 0.6 1.003 686 642
sd_participant__induced.statestress:Sexwoman:Drugplacebo 0.2 0.16 0.0 0.5 1.002 909 1888
sd_participant:stage__Intercept 0.1 0.06 0.0 0.2 1.004 621 1148
sd_session__Intercept 0.3 0.21 0.1 0.6 1.008 1041 1451
model: response.c ~ induced.state:Sex:Drug + (1 | session) + (1 | participant/stage) + (0 + induced.state:Sex:Drug | participant)
b_ = fix effect parameters
sd_ standard deviation of hierarchical parameters
ess = effective sample size


We calculate the same contrasts as above:

this time on from the posterior predictions:

fn = "brmsP1/stress_relief_contrasts.Rdata"
if (!file.exists(fn)) {
  postpre = unique(the_data$stage) %>% as.character() %>% sort()
  new_data = do.call(rbind,lapply(1:4, function(x) data.table(fit$data)[, session := x]))
  epred = posterior_epred(fit, newdata = new_data)
  obs.ratings = the_data[,.(r = mean(response)), by = .(response.c)][,r] %>% sort()
  epred.r = 
    apply(epred, 1, function(a) apply(a, 1, function(p) sum(p*obs.ratings)))
  
  epred.data = 
    new_data %>% 
    .[, .(induced.state, stage, Sex, participant, Drug)] %>% 
    .[, state.condition := ifelse(any(induced.state == "stress"),"stress","control"), 
      by = .(participant,Drug,induced.state)] %>% 
    .[, induced.state := NULL]
  
  pp = 
    cbind(epred.data, epred.r) %>% 
    melt(id.vars = names(epred.data), variable.name = "iter") %>% 
    .[, .(value = mean(value)), by = .(Sex,Drug,state.condition, iter)] %>% 
    .[, Sex := gsub("man","men",Sex)]
  
  
  # posterior distribution of contrasts, i.e.
  # main effects in men and women + their average and difference
  # main effects in men and women
  px = pp %>% 
    dcast(iter + Sex ~ Drug + state.condition) %>% 
    .[, stress_x_drug := 
        (oxycodone_stress - oxycodone_control) - 
        (placebo_stress - placebo_control)] %>% 
    .[, stress := 
        ((oxycodone_stress - oxycodone_control) + 
        (placebo_stress - placebo_control))/2] %>% 
    .[, .(iter, Sex, stress_x_drug, stress)]
  
  px = rbind(
    px,
    # average
    px[, .(stress_x_drug = mean(stress_x_drug),
           stress = mean(stress)),
       by = .(iter)][, Sex := "total"],
    # sex difference
    px %>% dcast(iter ~ Sex, value.var = c("stress","stress_x_drug")) %>% 
      .[, `:=`(stress_x_drug = stress_x_drug_men-stress_x_drug_women,
               stress = stress_men-stress_women)] %>% 
      .[, Sex := "Sex difference"] %>% 
      .[, .(iter, stress, stress_x_drug, Sex)]
  ) %>% 
    .[, Sex := factor(Sex, levels = c("men","women","Sex difference","total"))]
  save(px, pp, file = fn)
} else {
  load(fn)
}

Consistent with the simple stress effects seen in the raw data, we a credible stress effect is visible only in women.

p1 = pp %>% 
  dcast(iter + Sex + Drug ~ state.condition) %>% 
  .[, stress_effect := stress-control] %>% 
  ggplot(aes(x = Sex, y = stress_effect, color = Drug)) +
  stat_pointinterval(position = position_dodge(.2)) + 
  geom_hline(yintercept = 0, col = "red", lty = 2) + 
  scale_col_drug 

p2 = px %>% 
  ggplot(aes(x = stress_x_drug, fill = Sex)) + 
  geom_vline(xintercept = 0) + 
  geom_density(alpha = .5, color = NA) + 
  facet_wrap(~Sex) + 
  scale_col_sex.b

(p1 | p2) + 
  plot_layout(widths = c(2,3)) + 
  plot_annotation(tag_levels = "A")
A: Model-predicted effects of stress induction on 'stressed' ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Stress by Drug interaction in women, men, and sex difference (women-men)

Figure 7.2: A: Model-predicted effects of stress induction on ‘stressed’ ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Stress by Drug interaction in women, men, and sex difference (women-men)

We found no reliable effect of drug on stress. Oxycodone reduced stress induction induced increase of stress rating from pre state induction to post drug administration on average by NaN (CI=[NA, NA]; BFNA; P(O>P)NA). The interaction effect was -1 (CI=[-7, 6]; BF=0.7; P(O>P)=0.4) in women and -1 (CI=[-6, 5]; BF=0.6; P(O>P)=0.38) in men.

8 Buffering effects for reminders

Prepare the data:

my_data.Q_buffer = 
  my_data.Q %>% 
  .[stage %in% c(5,6) & 
      item_name %in% c(negaffect_items,
                       posaffect_items,
                       embselfcn_items,
                       "stressed")] %>% 
  .[, stress_reminder := ifelse(stage %in% c(5,7) & state.condition == "stress","pre","post")] %>% 
  .[, stress_reminder := factor(stress_reminder, levels = c("pre","post"))] %>% 
  .[, reminder := ifelse(stage %in% c(5,7),"pre","post")] %>% 
  .[, reminder := factor(reminder, levels = c("pre","post"))] %>% 
  .[, affect := ifelse(item_name %in% negaffect_items,
                       "negative",
                       ifelse(item_name %in% posaffect_items,
                              "positive","emb_selfc"))] %>% 
  .[item_name == "stressed", affect := "stressed"] %>% 
  coarsen(var = "response",
          bins = 20,
          range = c(0,100))

Here is a plot of the (nearly) raw data.

my_data.Q_buffer_stats = 
  my_data.Q_buffer %>% 
  .[, .(response = mean(response),
        se = sd(response)/sqrt(.N)),
    by = .(Drug, reminder, item_name, Sex, state.condition)] %>% 
  .[, `:=`(lower = response-2*se,
           upper = response+2*se)] %>% 
  .[, grp := paste(Drug,state.condition)]

my_data.Q_buffer_stats %>% 
  ggplot(aes(x = reminder, y = response, color = Drug, lty = state.condition, 
             shape = state.condition, group = grp,
             ymin = lower, ymax = upper)) +
  geom_point(position = pd.2) + 
  geom_line(aes(group = grp), position = pd.2) +
  geom_linerange(position = pd.2) + 
  facet_grid(item_name~Sex, scales = "free_y") +
  scale_col_drug
Mean ratings for stress, embarrassment and selfconciousness pre and post reminders by drug, state condition and sex.

Figure 8.1: Mean ratings for stress, embarrassment and selfconciousness pre and post reminders by drug, state condition and sex.

pdata = 
  my_data.Q_buffer %>% 
  # average the pre stages BL1, BL2, BL3
  .[, .(response = mean(response)), 
    by = .(participant, Drug, state.condition, stage, Sex,item_name )] %>% 
  # calculate state induction effect as post-pre difference
  dcast(participant  + Drug + state.condition + Sex + item_name ~ stage, value.var = "response") %>%
  .[, induction_effect := `6` - `5`] %>% 
  # calculate stress effect as difference of post-pre difference in stress & control conditions
  dcast(participant + Drug + Sex + item_name ~ state.condition, value.var = "induction_effect") %>% 
  .[, stress_effect := stress - control]

pdata.stats = 
  seWithin(pdata, "stress_effect",
           withinvars = c("Drug","item_name"),
           betweenvars = c("Sex"),
           idvar = "participant") %>% 
  data.table() %>% 
  .[, `:=`(lower=stress_effect-ci,
           upper=stress_effect+ci)]

pdata.stats %>% 
  ggplot(aes(x = Drug, y = stress_effect, fill = Drug)) + 
  geom_bar(stat = "identity", alpha = .5) + 
  geom_errorbar(aes(ymin = lower, ymax = upper), width = .25, linewidth = 1) + 
  facet_grid(item_name~Sex) + 
  geom_quasirandom(data = pdata, aes(color = Drug)) +
  scale_col_drug +
  ylim(-50,50)
Effect of stress reminders in the raw data measured after drug adminstration. Reminder effects are calculated as the change in stress ratings from pre to post stress induction and drug administration.

Figure 8.2: Effect of stress reminders in the raw data measured after drug adminstration. Reminder effects are calculated as the change in stress ratings from pre to post stress induction and drug administration.

pdata %>% 
  .[, stress_effect := stress - control] %>% 
  # calculate drug x stress interaction effect as difference of stress effects between drug conditions
  dcast(participant + item_name + Sex ~ Drug, value.var = "stress_effect") %>% 
  .[, stress_x_drug := oxycodone - placebo] %>% 
  ggplot(aes(x = stress_x_drug, fill=Sex)) + 
  geom_histogram(alpha = .5, position = "identity", bins = 20) + 
  scale_col_sex + 
  facet_wrap(~item_name)
Distribution of the interaction of stress effect (pre-post reminder) with drug condition (directly calculated from the data).

Figure 8.3: Distribution of the interaction of stress effect (pre-post reminder) with drug condition (directly calculated from the data).

8.1 Analysis approach

The analysis approach here is that we get posterior distributions for the ratings in the 8 cells obtained by crossing:

  • stress vs control
  • pre vs post reminder
  • drug vs placebo
  • males vs females

Then we can calculate the desired contrasts by first calculating the pre to post change for 4 conditions from stress x drug. - pre_post_placebo_control = Drugplacebo:Stresscontrol:reminderpre - Drugplacebo:Stresscontrol:reminderpost - pre_post_placebo_stress = Drugplacebo:Stressstress:reminderpre - Drugplacebo:Stressstress:reminderpost - pre_post_placebo_control = Drugplacebo:Stresscontrol:reminderpre - Drugplacebo:Stresscontrol:reminderpost - pre_post_placebo_stress = Drugplacebo:Stressstress:reminderpre - Drugplacebo:Stressstress:reminderpost

With those, we can calculate the stress reminder effects in the drug and placebo condition

  • stress_effect_placebo = pre_post_placebo_stress - pre_post_placebo_control
  • stress_effect_placebo = pre_post_placebo_stress - pre_post_placebo_control

And finally we can calculate if the stress-reminder effect is smaller in the the placebo compared to the placebo condition

  • stress_x_drug = stress_effect_placebo - stress_effect_placebo

This description does not include sex effects, but these effects can be calculates separately by sex.

Note that given this specification of the contrasts, a buffer effect where post-reminder ratings of negative items like stress are lower under oxycodone compared to placebo will result in negative values for the stress_x_drug interaction effect.

Here is the analysis model:

analysis_formula = 
  response.c ~  Sex:Drug:Stress:reminder + (1|session) + (1 + Drug:Stress:reminder | participant) 

The analysis for manipulation checks showed that that the stress effects of the first reminder were clearer, especially for feeling stressed. Therefore, the analysis of buffer effects focuses on the first reminder.

8.2 Buffering of stress

We focus on buffering of stress for the first reminder.

f = "brmsP1/buffer/stressed.Rdata"
the_data = my_data.Q_buffer[item_name == "stressed"]
if (file.exists(f)) {
  load(f)
} else {
  fit = fit_ordinal_model(
    the_data, analysis_formula, f, family = acat())
}
Click to show sample diagnostics
Table 8.1: Summary of model paramters
parameter mean sd q5 q95 rhat ess_bulk ess_tail
Threshold[1] -2.0 0.97 -3.6 -0.4 1.014 394 409
Threshold[2] -1.7 0.88 -3.2 -0.3 1.013 452 521
Threshold[3] -0.8 0.84 -2.2 0.6 1.013 500 755
Threshold[4] -0.6 0.84 -2.0 0.8 1.012 566 801
Threshold[5] -0.3 0.84 -1.7 1.1 1.011 624 938
Threshold[6] 0.0 0.84 -1.4 1.4 1.012 646 1000
Threshold[7] 0.5 0.86 -0.9 1.9 1.010 770 1494
Threshold[8] 0.3 0.87 -1.2 1.7 1.010 726 1568
Threshold[9] 0.8 0.86 -0.7 2.2 1.008 828 1672
Threshold[10] 1.6 0.93 0.0 3.1 1.008 987 1638
Threshold[11] 1.4 0.99 -0.2 3.1 1.006 1091 2064
Threshold[12] 3.0 1.22 1.0 5.0 1.003 1511 2240
Threshold[13] -0.5 1.17 -2.4 1.4 1.006 1385 2313
Threshold[14] 3.4 1.08 1.6 5.2 1.006 1327 2325
Threshold[15] 1.9 1.16 -0.1 3.8 1.004 1552 2328
Threshold[16] 3.1 1.17 1.2 5.1 1.004 1440 2207
Threshold[17] 3.5 1.27 1.4 5.6 1.004 1724 2671
b_Sexman:Drugoxycodone:Stresscontrol:reminderpre -1.1 0.72 -2.3 0.1 1.007 894 1487
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpre -0.5 0.68 -1.7 0.6 1.008 883 1976
b_Sexman:Drugplacebo:Stresscontrol:reminderpre -0.8 0.63 -1.9 0.2 1.007 774 1486
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpre -0.2 0.62 -1.2 0.8 1.007 869 1573
b_Sexman:Drugoxycodone:Stressstress:reminderpre -0.2 0.66 -1.3 0.9 1.004 938 1380
b_Sexwoman:Drugoxycodone:Stressstress:reminderpre 0.7 0.65 -0.4 1.7 1.007 905 1500
b_Sexman:Drugplacebo:Stressstress:reminderpre 0.2 0.60 -0.8 1.2 1.004 895 1606
b_Sexwoman:Drugplacebo:Stressstress:reminderpre 0.8 0.60 -0.2 1.8 1.005 857 1584
b_Sexman:Drugoxycodone:Stresscontrol:reminderpost -1.0 0.77 -2.3 0.3 1.005 1162 1845
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpost -0.6 0.75 -1.8 0.6 1.005 1081 2023
b_Sexman:Drugplacebo:Stresscontrol:reminderpost -1.2 0.67 -2.3 -0.1 1.007 834 1387
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpost -0.6 0.67 -1.8 0.4 1.005 943 1786
b_Sexman:Drugoxycodone:Stressstress:reminderpost 0.8 0.66 -0.3 1.9 1.007 882 1465
b_Sexwoman:Drugoxycodone:Stressstress:reminderpost 1.6 0.67 0.5 2.7 1.007 848 1549
b_Sexman:Drugplacebo:Stressstress:reminderpost 0.1 0.65 -1.0 1.2 1.004 973 1655
b_Sexwoman:Drugplacebo:Stressstress:reminderpost 1.9 0.66 0.8 3.0 1.007 831 1539
sd_participant__Intercept 1.5 0.32 1.0 2.0 1.003 697 1364
sd_participant__Drugoxycodone:Stresscontrol:reminderpre 1.9 0.48 1.2 2.8 1.003 606 1339
sd_participant__Drugplacebo:Stresscontrol:reminderpre 0.7 0.42 0.1 1.5 1.005 435 1106
sd_participant__Drugoxycodone:Stressstress:reminderpre 1.6 0.43 1.0 2.4 1.006 591 1150
sd_participant__Drugplacebo:Stressstress:reminderpre 1.0 0.39 0.4 1.7 1.013 410 566
sd_participant__Drugoxycodone:Stresscontrol:reminderpost 2.2 0.50 1.5 3.1 1.010 527 877
sd_participant__Drugplacebo:Stresscontrol:reminderpost 1.1 0.47 0.4 2.0 1.013 389 440
sd_participant__Drugoxycodone:Stressstress:reminderpost 2.0 0.46 1.3 2.8 1.003 435 720
sd_participant__Drugplacebo:Stressstress:reminderpost 1.5 0.44 0.8 2.3 1.009 380 535
sd_session__Intercept 2.6 0.93 1.2 4.2 1.003 632 405
model: response.c ~ Sex:Drug:Stress:reminder + (1 | session) + (1 + Drug:Stress:reminder | participant)
b_ = fix effect parameters
sd_ standard deviation of hierarchical parameters
ess = effective sample size


p = plot_buffer_effect(fit,the_data,"stressed")
p$ppc
Posterior predictive check of pre-post differences. Points are differences directly calculate from the data, intervals show differences calculate from model parameters.

Figure 8.4: Posterior predictive check of pre-post differences. Points are differences directly calculate from the data, intervals show differences calculate from model parameters.

The following figure shows stress-drug interactions (negative values indicate a buffering effect). The data support a hypothesis of no effect of oxycodone on stress response if the credible interval falls in a region of practical equivalence of +/- 5 points.

((p$stress.plot + 
    ylab("stress effect on 'stressed' ratings.")) | p$buffer.plot) + 
  plot_annotation(tag_levels = "A") + 
  plot_layout(widths = c(2,3))
A: Model-predicted effects of the first stress reminder on 'stressed' ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Posterior distribution of the interaction of stress effect (pre-post reminder) on 'stressed' ratings with drug condition. Dashed vertical lines enclode a region of practical equivalence.

Figure 8.5: A: Model-predicted effects of the first stress reminder on ‘stressed’ ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Posterior distribution of the interaction of stress effect (pre-post reminder) on ‘stressed’ ratings with drug condition. Dashed vertical lines enclode a region of practical equivalence.

p$buffer.stats %>% 
  kable(caption = "Buffering effect of oxycodone for 'stressed' ratings after reminders.") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic()
Table 8.2: Buffering effect of oxycodone for ‘stressed’ ratings after reminders.
Sex stress stress in placebo stress in oxycodone buffer
men 4 [2, 7]; >.999 4 [0, 7]; 0.9775 5 [1, 9]; 0.997 2 [-4, 7]; 0.72725
women 12 [9, 15]; >.999 16 [12, 20]; >.999 8 [3, 12]; >.999 -8 [-14, -2]; 0.0035
total 8 [6, 10]; >.999 10 [7, 13]; >.999 6 [4, 9]; >.999 -3 [-7, 1]; 0.05
Sex difference 7 [3, 11]; >.999 12 [6, 18]; >.999 2 [-3, 8]; 0.7955 -10 [-18, -2]; 0.0095

8.3 Buffering of embarrassment

f = "brmsP1/buffer/emb.Rdata"
the_data = my_data.Q_buffer[item_name == "embarassed"]
if (file.exists(f)) {
  load(f)
} else {
  fit = fit_ordinal_model(
    the_data, analysis_formula, f, family = acat())
}
Click to show sample diagnostics
Table 8.3: Summary of model paramters
parameter mean sd q5 q95 rhat ess_bulk ess_tail
Threshold[1] -2.0 0.87 -3.4 -0.5 1.003 920 1447
Threshold[2] -1.5 0.80 -2.8 -0.2 1.002 990 1480
Threshold[3] -1.2 0.78 -2.5 0.1 1.001 1067 1540
Threshold[4] -0.5 0.79 -1.8 0.8 1.002 1113 1589
Threshold[5] -0.2 0.81 -1.5 1.2 1.002 1237 1720
Threshold[6] 0.2 0.84 -1.2 1.6 1.001 1265 1890
Threshold[7] 0.5 0.90 -0.9 2.0 1.002 1628 2181
Threshold[8] 0.5 0.94 -1.0 2.0 1.001 1682 1907
Threshold[9] -0.5 0.91 -2.1 1.0 1.001 1517 2289
Threshold[10] 2.0 0.95 0.4 3.5 1.000 1745 2642
Threshold[11] 0.8 1.06 -0.9 2.6 1.000 1966 2063
Threshold[12] 0.3 1.00 -1.3 2.0 1.001 1901 2965
Threshold[13] 1.3 0.90 -0.1 2.8 1.000 1521 2469
Threshold[14] 1.3 0.88 -0.2 2.7 1.001 1576 2248
Threshold[15] 4.3 1.20 2.4 6.3 1.002 2299 2466
Threshold[16] 1.4 1.28 -0.6 3.5 1.000 2606 2651
Threshold[17] 2.2 1.15 0.3 4.1 1.002 2114 2831
Threshold[18] 2.5 1.08 0.8 4.3 1.000 1859 2689
b_Sexman:Drugoxycodone:Stresscontrol:reminderpre -1.1 0.73 -2.4 0.0 1.001 1552 2509
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpre -0.4 0.70 -1.6 0.8 1.001 1624 2183
b_Sexman:Drugplacebo:Stresscontrol:reminderpre -0.2 0.71 -1.4 1.0 1.001 1723 2268
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpre -1.3 0.81 -2.7 0.0 1.002 1760 1997
b_Sexman:Drugoxycodone:Stressstress:reminderpre -0.3 0.73 -1.5 0.9 1.002 1590 2467
b_Sexwoman:Drugoxycodone:Stressstress:reminderpre 0.3 0.75 -1.0 1.5 1.001 1533 1957
b_Sexman:Drugplacebo:Stressstress:reminderpre -0.4 0.73 -1.6 0.8 1.001 1805 2477
b_Sexwoman:Drugplacebo:Stressstress:reminderpre 1.3 0.68 0.1 2.4 1.001 1545 2003
b_Sexman:Drugoxycodone:Stresscontrol:reminderpost -0.7 0.71 -1.9 0.4 1.002 1521 2407
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpost -1.6 0.78 -2.9 -0.3 1.003 1421 2283
b_Sexman:Drugplacebo:Stresscontrol:reminderpost -1.0 0.67 -2.1 0.1 1.001 1572 2120
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpost -1.6 0.75 -2.8 -0.4 1.001 1528 2385
b_Sexman:Drugoxycodone:Stressstress:reminderpost 0.9 0.71 -0.2 2.1 1.001 1640 2468
b_Sexwoman:Drugoxycodone:Stressstress:reminderpost 1.6 0.72 0.4 2.8 1.001 1555 2074
b_Sexman:Drugplacebo:Stressstress:reminderpost 0.9 0.67 -0.1 2.0 1.002 1582 2393
b_Sexwoman:Drugplacebo:Stressstress:reminderpost 3.5 0.75 2.3 4.8 1.002 1254 2026
sd_participant__Intercept 1.9 0.39 1.4 2.6 1.002 1303 1937
sd_participant__Drugoxycodone:Stresscontrol:reminderpre 1.3 0.48 0.5 2.1 1.008 1050 965
sd_participant__Drugplacebo:Stresscontrol:reminderpre 1.7 0.45 1.0 2.5 1.002 1397 2582
sd_participant__Drugoxycodone:Stressstress:reminderpre 2.1 0.50 1.3 2.9 1.002 1306 2050
sd_participant__Drugplacebo:Stressstress:reminderpre 1.8 0.45 1.1 2.5 1.010 783 1111
sd_participant__Drugoxycodone:Stresscontrol:reminderpost 1.0 0.50 0.2 1.8 1.014 662 1056
sd_participant__Drugplacebo:Stresscontrol:reminderpost 0.6 0.44 0.1 1.5 1.007 816 1699
sd_participant__Drugoxycodone:Stressstress:reminderpost 2.1 0.50 1.4 3.0 1.002 1102 1792
sd_participant__Drugplacebo:Stressstress:reminderpost 1.8 0.44 1.1 2.6 1.015 745 1178
sd_session__Intercept 3.0 0.93 1.7 4.7 1.003 1435 1728
model: response.c ~ Sex:Drug:Stress:reminder + (1 | session) + (1 + Drug:Stress:reminder | participant)
b_ = fix effect parameters
sd_ standard deviation of hierarchical parameters
ess = effective sample size


p = plot_buffer_effect(fit,the_data,"embarassed")
p$ppc
Posterior predictive check of pre-post differences. Points are differences directly calculate from the data, intervals show differences calculate from model parameters.

Figure 8.6: Posterior predictive check of pre-post differences. Points are differences directly calculate from the data, intervals show differences calculate from model parameters.

((p$stress.plot + 
    ylab("stress effect on 'embarassed' ratings.")) | p$buffer.plot) + 
  plot_annotation(tag_levels = "A") + 
  plot_layout(widths = c(2,3))
A: Model-predicted effects of the first stress reminder on 'embarassed' ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Posterior distribution of the interaction of stress effect (pre-post reminder) on 'embarassed' ratings with drug condition. Dashed vertical lines enclode a region of practical equivalence.

Figure 8.7: A: Model-predicted effects of the first stress reminder on ‘embarassed’ ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Posterior distribution of the interaction of stress effect (pre-post reminder) on ‘embarassed’ ratings with drug condition. Dashed vertical lines enclode a region of practical equivalence.

p$buffer.stats %>% 
  kable(caption = "Buffering effect of oxycodone for 'embarassed' ratings after reminders.") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic()
Table 8.4: Buffering effect of oxycodone for ‘embarassed’ ratings after reminders.
Sex stress stress in placebo stress in oxycodone buffer
men 9 [7, 12]; >.999 11 [7, 14]; >.999 7 [4, 11]; >.999 -3 [-8, 2]; 0.09525
women 22 [20, 25]; >.999 30 [25, 34]; >.999 15 [11, 19]; >.999 -15 [-20, -9]; <.001
total 16 [14, 18]; >.999 20 [17, 23]; >.999 11 [9, 14]; >.999 -9 [-13, -5]; <.001
Sex difference 13 [10, 17]; >.999 19 [13, 25]; >.999 8 [3, 13]; 0.9985 -11 [-19, -3]; 0.00225

8.4 Buffering of selfconciousness

f = "brmsP1/buffer/selfcn.Rdata"
the_data = my_data.Q_buffer[item_name == "selfconscious"]
if (file.exists(f)) {
  load(f)
} else {
  fit = fit_ordinal_model(
    the_data, analysis_formula, f, family = acat())
}
Click to show sample diagnostics
Table 8.5: Summary of model paramters
parameter mean sd q5 q95 rhat ess_bulk ess_tail
Threshold[1] -2.3 0.78 -3.5 -1.0 1.012 500 1004
Threshold[2] -1.6 0.80 -3.0 -0.3 1.008 672 1286
Threshold[3] -3.2 0.80 -4.5 -1.9 1.007 703 1569
Threshold[4] -1.8 0.74 -3.0 -0.6 1.011 648 1189
Threshold[5] -1.8 0.73 -3.0 -0.6 1.009 665 1429
Threshold[6] -2.2 0.70 -3.4 -1.1 1.010 609 1349
Threshold[7] -1.3 0.66 -2.4 -0.3 1.010 603 1289
Threshold[8] -0.3 0.71 -1.5 0.9 1.007 707 1567
Threshold[9] -2.5 0.67 -3.7 -1.5 1.008 655 1381
Threshold[10] 0.1 0.61 -0.9 1.0 1.008 581 1214
Threshold[11] -0.1 0.63 -1.2 0.9 1.009 585 1271
Threshold[12] -0.1 0.61 -1.2 0.8 1.006 578 1219
Threshold[13] 0.8 0.63 -0.2 1.8 1.006 571 1200
Threshold[14] 0.8 0.64 -0.3 1.8 1.006 542 1124
Threshold[15] 1.8 0.66 0.7 2.9 1.006 555 1182
Threshold[16] 1.2 0.69 0.1 2.4 1.002 522 1044
Threshold[17] 2.3 0.70 1.2 3.5 1.004 491 964
Threshold[18] 3.5 0.76 2.2 4.7 1.002 502 1011
Threshold[19] 2.7 0.83 1.3 4.0 1.003 450 909
b_Sexman:Drugoxycodone:Stresscontrol:reminderpre 0.0 0.64 -1.0 1.1 1.009 453 867
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpre -0.4 0.61 -1.4 0.6 1.006 523 987
b_Sexman:Drugplacebo:Stresscontrol:reminderpre 0.0 0.62 -1.0 1.0 1.008 449 918
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpre -0.2 0.59 -1.1 0.8 1.007 505 1061
b_Sexman:Drugoxycodone:Stressstress:reminderpre 0.2 0.63 -0.9 1.2 1.006 477 1138
b_Sexwoman:Drugoxycodone:Stressstress:reminderpre -0.5 0.60 -1.5 0.5 1.005 575 1181
b_Sexman:Drugplacebo:Stressstress:reminderpre 0.4 0.63 -0.6 1.4 1.006 481 836
b_Sexwoman:Drugplacebo:Stressstress:reminderpre 0.1 0.60 -0.8 1.1 1.005 510 969
b_Sexman:Drugoxycodone:Stresscontrol:reminderpost 0.4 0.62 -0.7 1.4 1.007 427 867
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpost -0.5 0.60 -1.5 0.4 1.006 510 1034
b_Sexman:Drugplacebo:Stresscontrol:reminderpost -0.2 0.65 -1.3 0.9 1.007 437 904
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpost -0.8 0.63 -1.8 0.3 1.006 544 1040
b_Sexman:Drugoxycodone:Stressstress:reminderpost 0.3 0.61 -0.7 1.3 1.006 495 926
b_Sexwoman:Drugoxycodone:Stressstress:reminderpost -0.2 0.58 -1.1 0.8 1.006 550 1158
b_Sexman:Drugplacebo:Stressstress:reminderpost 0.5 0.64 -0.6 1.5 1.006 512 1014
b_Sexwoman:Drugplacebo:Stressstress:reminderpost 1.0 0.63 0.0 2.0 1.004 559 954
sd_participant__Intercept 2.4 0.41 1.8 3.1 1.010 406 923
sd_participant__Drugoxycodone:Stresscontrol:reminderpre 1.4 0.36 0.9 2.0 1.015 296 341
sd_participant__Drugplacebo:Stresscontrol:reminderpre 0.7 0.38 0.1 1.4 1.027 154 287
sd_participant__Drugoxycodone:Stressstress:reminderpre 1.5 0.33 1.0 2.1 1.009 364 871
sd_participant__Drugplacebo:Stressstress:reminderpre 0.9 0.34 0.4 1.5 1.010 150 355
sd_participant__Drugoxycodone:Stresscontrol:reminderpost 1.1 0.33 0.6 1.7 1.014 260 235
sd_participant__Drugplacebo:Stresscontrol:reminderpost 1.3 0.35 0.7 1.9 1.007 269 944
sd_participant__Drugoxycodone:Stressstress:reminderpost 1.2 0.30 0.8 1.8 1.018 368 772
sd_participant__Drugplacebo:Stressstress:reminderpost 1.6 0.36 1.1 2.2 1.004 259 861
sd_session__Intercept 0.1 0.16 0.0 0.4 1.001 1390 1909
model: response.c ~ Sex:Drug:Stress:reminder + (1 | session) + (1 + Drug:Stress:reminder | participant)
b_ = fix effect parameters
sd_ standard deviation of hierarchical parameters
ess = effective sample size


p = plot_buffer_effect(fit,the_data,"selfconscious")
p$ppc
Posterior predictive check of pre-post differences. Points are differences directly calculate from the data, intervals show differences calculate from model parameters.

Figure 8.8: Posterior predictive check of pre-post differences. Points are differences directly calculate from the data, intervals show differences calculate from model parameters.

((p$stress.plot + 
    ylab("stress effect on 'selfconscious' ratings.")) | p$buffer.plot) + 
  plot_annotation(tag_levels = "A") + 
  plot_layout(widths = c(2,3))
A: Model-predicted effects of the first stress reminder on 'selfconscious' ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Posterior distribution of the interaction of stress effect (pre-post reminder) on 'selfconscious' ratings with drug condition. Dashed vertical lines enclode a region of practical equivalence.

Figure 8.9: A: Model-predicted effects of the first stress reminder on ‘selfconscious’ ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Posterior distribution of the interaction of stress effect (pre-post reminder) on ‘selfconscious’ ratings with drug condition. Dashed vertical lines enclode a region of practical equivalence.

p$buffer.stats %>% 
  kable(caption = "Buffering effect of oxycodone for 'selfconscious' ratings after reminders.") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic()
Table 8.6: Buffering effect of oxycodone for ‘selfconscious’ ratings after reminders.
Sex stress stress in placebo stress in oxycodone buffer
men 1 [-3, 5]; 0.73425 4 [-1, 10]; 0.94775 -2 [-7, 3]; 0.244 -6 [-13, 1]; 0.051
women 12 [8, 16]; >.999 18 [12, 24]; >.999 6 [0, 12]; 0.97675 -12 [-20, -4]; 0.0025
total 7 [4, 9]; >.999 11 [7, 15]; >.999 2 [-2, 6]; 0.85225 -9 [-14, -4]; <.001
Sex difference 11 [5, 17]; >.999 14 [6, 22]; >.999 8 [-0, 16]; 0.97175 -6 [-17, 5]; 0.14725

8.5 Buffering of stress, embarrassement & selfconciousness

Given that the data patterns for these three items are expected to be similar, we can obtain more stable estimates of a buffer effect by modeling all items jointly.

f = "brmsP1/buffer/embselfstress.Rdata"
the_data = my_data.Q_buffer[item_name %in% c(embselfcn_items,"stressed")]
analysis_formula = 
  response.c ~  Sex:Drug:Stress:reminder + (1|session) + (Drug:Stress:reminder| participant) + (1|item_name)
if (file.exists(f)) {
  load(f)
} else {
  fit = fit_ordinal_model(
    the_data, analysis_formula, f, family = acat())
}
Click to show sample diagnostics
Table 8.7: Summary of model paramters
parameter mean sd q5 q95 rhat ess_bulk ess_tail
Threshold[1] 1.1 0.54 0.3 2.0 1.007 631 753
Threshold[2] 0.3 0.55 -0.6 1.2 1.010 648 918
Threshold[3] 0.1 0.56 -0.8 1.0 1.005 660 956
Threshold[4] 0.3 0.56 -0.7 1.2 1.008 674 1084
Threshold[5] 0.2 0.58 -0.8 1.1 1.006 711 1253
Threshold[6] -0.1 0.58 -1.0 0.9 1.005 697 1070
Threshold[7] 0.2 0.58 -0.7 1.2 1.006 724 1178
Threshold[8] 0.5 0.59 -0.5 1.5 1.004 768 1340
Threshold[9] -0.9 0.58 -1.9 0.0 1.007 707 1186
Threshold[10] 0.9 0.57 -0.1 1.8 1.006 698 1089
Threshold[11] 0.3 0.58 -0.6 1.3 1.006 747 898
Threshold[12] 0.0 0.59 -0.9 1.0 1.007 701 1265
Threshold[13] 0.4 0.58 -0.6 1.3 1.004 699 888
Threshold[14] 0.3 0.58 -0.6 1.3 1.006 721 1061
Threshold[15] 1.0 0.60 0.0 2.0 1.007 726 1379
Threshold[16] 0.0 0.59 -1.0 1.0 1.006 753 1050
Threshold[17] 0.5 0.58 -0.5 1.4 1.005 709 998
Threshold[18] 1.3 0.61 0.3 2.3 1.003 757 1502
Threshold[19] -0.6 0.61 -1.6 0.4 1.006 764 1183
b_Sexman:Drugoxycodone:Stresscontrol:reminderpre -0.1 0.48 -0.8 0.7 1.010 615 1092
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpre -0.1 0.47 -0.8 0.8 1.010 617 987
b_Sexman:Drugplacebo:Stresscontrol:reminderpre 0.0 0.48 -0.8 0.8 1.010 615 1067
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpre -0.1 0.47 -0.8 0.7 1.010 616 977
b_Sexman:Drugoxycodone:Stressstress:reminderpre 0.0 0.48 -0.8 0.8 1.010 614 1128
b_Sexwoman:Drugoxycodone:Stressstress:reminderpre 0.0 0.47 -0.7 0.8 1.010 616 979
b_Sexman:Drugplacebo:Stressstress:reminderpre 0.0 0.48 -0.7 0.8 1.009 613 1077
b_Sexwoman:Drugplacebo:Stressstress:reminderpre 0.1 0.47 -0.7 0.9 1.010 615 1045
b_Sexman:Drugoxycodone:Stresscontrol:reminderpost 0.0 0.47 -0.8 0.8 1.011 619 1096
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpost -0.1 0.47 -0.8 0.7 1.011 623 1038
b_Sexman:Drugplacebo:Stresscontrol:reminderpost -0.1 0.48 -0.9 0.7 1.011 616 1095
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpost -0.1 0.47 -0.9 0.7 1.011 616 985
b_Sexman:Drugoxycodone:Stressstress:reminderpost 0.1 0.48 -0.7 0.9 1.009 621 1128
b_Sexwoman:Drugoxycodone:Stressstress:reminderpost 0.1 0.47 -0.6 0.9 1.011 616 1012
b_Sexman:Drugplacebo:Stressstress:reminderpost 0.1 0.48 -0.7 0.9 1.010 622 1093
b_Sexwoman:Drugplacebo:Stressstress:reminderpost 0.3 0.47 -0.4 1.1 1.011 620 1042
sd_item_name__Intercept 0.5 0.30 0.2 1.2 1.003 1446 1774
sd_participant__Intercept 0.2 0.03 0.2 0.3 1.006 1391 1652
sd_participant__Drugoxycodone:Stresscontrol:reminderpre 0.1 0.04 0.0 0.1 1.004 1010 1664
sd_participant__Drugplacebo:Stresscontrol:reminderpre 0.0 0.03 0.0 0.1 1.001 2222 2015
sd_participant__Drugoxycodone:Stressstress:reminderpre 0.0 0.03 0.0 0.1 1.001 1470 1557
sd_participant__Drugplacebo:Stressstress:reminderpre 0.0 0.03 0.0 0.1 1.000 1458 1865
sd_participant__Drugoxycodone:Stresscontrol:reminderpost 0.1 0.04 0.0 0.1 1.003 1262 1970
sd_participant__Drugplacebo:Stresscontrol:reminderpost 0.0 0.03 0.0 0.1 1.001 2456 2618
sd_participant__Drugoxycodone:Stressstress:reminderpost 0.2 0.04 0.1 0.2 1.001 1644 1711
sd_participant__Drugplacebo:Stressstress:reminderpost 0.2 0.04 0.1 0.2 1.002 1550 1778
sd_session__Intercept 0.0 0.04 0.0 0.1 1.000 1493 1782
model: response.c ~ Sex:Drug:Stress:reminder + (1 | session) + (Drug:Stress:reminder | participant) + (1 | item_name)
b_ = fix effect parameters
sd_ standard deviation of hierarchical parameters
ess = effective sample size


p = plot_buffer_effect(fit,the_data,"stress_emb_selfc")
((p$stress.plot + 
    ylab("stress effect on 'stressed', 'emb.' & 'selfc.' ratings.")) | p$buffer.plot) + 
  plot_annotation(tag_levels = "A") + 
  plot_layout(widths = c(2,3))
A: Model-predicted effects of the first stress reminder on 'stressed', 'selfconscious' and 'embarrased' ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Posterior distribution of the interaction of stress effect (pre-post reminder) ratings with drug condition. Dashed vertical lines enclode a region of practical equivalence.

Figure 8.10: A: Model-predicted effects of the first stress reminder on ‘stressed’, ‘selfconscious’ and ‘embarrased’ ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Posterior distribution of the interaction of stress effect (pre-post reminder) ratings with drug condition. Dashed vertical lines enclode a region of practical equivalence.

p$buffer.stats %>% 
  kable(caption = "Buffering effect of oxycodone for 'stressed', 'emb.' & 'selfc.' ratings after reminders.") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic()
Table 8.8: Buffering effect of oxycodone for ‘stressed’, ‘emb.’ & ‘selfc.’ ratings after reminders.
Sex stress stress in placebo stress in oxycodone buffer
men 5 [0, 10]; 0.97975 6 [-1, 13]; 0.9675 3 [-4, 10]; 0.83775 -3 [-13, 7]; 0.28525
women 15 [11, 20]; >.999 21 [14, 28]; >.999 10 [3, 16]; 0.9965 -12 [-22, -2]; 0.00925
total 10 [7, 13]; >.999 14 [9, 18]; >.999 6 [2, 11]; 0.9965 -7 [-14, -0]; 0.023
Sex difference 11 [4, 17]; 0.99875 15 [5, 25]; 0.998 6 [-4, 16]; 0.88525 -9 [-23, 5]; 0.098

8.6 Buffering of negative mood

f = "brmsP1/buffer/negaffect.Rdata"
the_data = my_data.Q_buffer[item_name %in% negaffect_items]
if (file.exists(f)) {
  load(f)
} else {
  fit = fit_ordinal_model(
    the_data, analysis_formula, f, family = acat())
}
Click to show sample diagnostics
Table 8.9: Summary of model paramters
parameter mean sd q5 q95 rhat ess_bulk ess_tail
Threshold[1] 1.0 0.60 0.0 2.0 1.002 851 1757
Threshold[2] 0.6 0.61 -0.4 1.6 1.001 864 1748
Threshold[3] 0.9 0.62 -0.1 1.9 1.001 874 1954
Threshold[4] 0.6 0.63 -0.4 1.7 1.002 929 1816
Threshold[5] 0.9 0.64 -0.2 1.9 1.001 922 1916
Threshold[6] 0.8 0.65 -0.3 1.8 1.002 966 1896
Threshold[7] 0.6 0.65 -0.5 1.6 1.002 996 1984
Threshold[8] 1.3 0.67 0.2 2.4 1.001 1052 2125
Threshold[9] 0.4 0.68 -0.8 1.5 1.001 1048 2193
Threshold[10] 1.0 0.67 -0.1 2.1 1.001 1018 2233
Threshold[11] 1.2 0.69 0.0 2.3 1.001 1097 2187
Threshold[12] 0.7 0.69 -0.4 1.9 1.001 1118 2420
Threshold[13] 1.5 0.71 0.3 2.6 1.002 1161 2360
Threshold[14] 1.2 0.75 0.0 2.5 1.001 1277 2810
Threshold[15] 2.2 0.88 0.8 3.7 1.001 1609 3250
Threshold[16] 0.1 0.87 -1.4 1.5 1.001 1632 3193
Threshold[17] 1.9 0.86 0.5 3.3 1.001 1620 3411
Threshold[18] 1.9 1.02 0.3 3.6 1.002 1860 3197
Threshold[19] -0.1 0.95 -1.7 1.4 1.001 2040 4168
b_Sexman:Drugoxycodone:Stresscontrol:reminderpre -0.4 0.54 -1.2 0.5 1.003 892 2056
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpre 0.1 0.52 -0.7 1.0 1.003 878 1978
b_Sexman:Drugplacebo:Stresscontrol:reminderpre -0.4 0.53 -1.2 0.5 1.004 867 1914
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpre -0.2 0.53 -1.0 0.7 1.004 845 1877
b_Sexman:Drugoxycodone:Stressstress:reminderpre -0.1 0.52 -1.0 0.7 1.003 861 1862
b_Sexwoman:Drugoxycodone:Stressstress:reminderpre 0.4 0.52 -0.4 1.3 1.003 842 1870
b_Sexman:Drugplacebo:Stressstress:reminderpre 0.0 0.52 -0.9 0.8 1.003 866 2070
b_Sexwoman:Drugplacebo:Stressstress:reminderpre 0.2 0.52 -0.6 1.1 1.004 845 1884
b_Sexman:Drugoxycodone:Stresscontrol:reminderpost -0.5 0.55 -1.4 0.4 1.003 941 2108
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpost 0.0 0.54 -0.9 0.9 1.003 909 2029
b_Sexman:Drugplacebo:Stresscontrol:reminderpost -0.6 0.54 -1.4 0.3 1.003 884 1853
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpost -0.3 0.53 -1.1 0.6 1.004 850 1883
b_Sexman:Drugoxycodone:Stressstress:reminderpost 0.0 0.52 -0.8 0.9 1.004 841 1906
b_Sexwoman:Drugoxycodone:Stressstress:reminderpost 0.6 0.52 -0.2 1.5 1.003 849 1748
b_Sexman:Drugplacebo:Stressstress:reminderpost 0.2 0.52 -0.6 1.1 1.003 848 1913
b_Sexwoman:Drugplacebo:Stressstress:reminderpost 0.6 0.52 -0.2 1.5 1.004 840 1911
sd_item_name__Intercept 0.6 0.38 0.2 1.3 1.003 1492 2131
sd_participant__Intercept 0.7 0.11 0.5 0.8 1.001 2427 4003
sd_participant__Drugoxycodone:Stresscontrol:reminderpre 0.6 0.14 0.3 0.8 1.005 1049 1800
sd_participant__Drugplacebo:Stresscontrol:reminderpre 0.3 0.19 0.0 0.7 1.005 403 1241
sd_participant__Drugoxycodone:Stressstress:reminderpre 0.4 0.12 0.2 0.6 1.003 526 648
sd_participant__Drugplacebo:Stressstress:reminderpre 0.4 0.14 0.2 0.6 1.005 659 753
sd_participant__Drugoxycodone:Stresscontrol:reminderpost 0.7 0.17 0.4 1.0 1.004 1263 2195
sd_participant__Drugplacebo:Stresscontrol:reminderpost 0.4 0.22 0.1 0.8 1.003 434 1152
sd_participant__Drugoxycodone:Stressstress:reminderpost 0.4 0.11 0.2 0.6 1.003 575 1427
sd_participant__Drugplacebo:Stressstress:reminderpost 0.4 0.13 0.2 0.6 1.004 522 662
sd_session__Intercept 0.2 0.23 0.1 0.7 1.001 2363 2642
model: response.c ~ Sex:Drug:Stress:reminder + (1 | session) + (Drug:Stress:reminder | participant) + (1 | item_name)
b_ = fix effect parameters
sd_ standard deviation of hierarchical parameters
ess = effective sample size


p = plot_buffer_effect(fit,the_data,"negaffect")
((p$stress.plot + 
    ylab("stress effect on negative mood ratings")) | p$buffer.plot) + 
  plot_annotation(tag_levels = "A") + 
  plot_layout(widths = c(2,3))
A: Model-predicted effects of the first stress reminder on 'negative mood', 'selfconscious' and 'embarrased' ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Posterior distribution of the interaction of stress effect (pre-post reminder) ratings with drug condition. Dashed vertical lines enclode a region of practical equivalence.

Figure 8.11: A: Model-predicted effects of the first stress reminder on ‘negative mood’, ‘selfconscious’ and ‘embarrased’ ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Posterior distribution of the interaction of stress effect (pre-post reminder) ratings with drug condition. Dashed vertical lines enclode a region of practical equivalence.

p$buffer.stats %>% 
  kable(caption = "Buffering effect of oxycodone for negative mood ratings after reminders.") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic()
Table 8.10: Buffering effect of oxycodone for negative mood ratings after reminders.
Sex stress stress in placebo stress in oxycodone buffer
men 3 [1, 6]; >.999 5 [2, 8]; >.999 2 [-1, 5]; 0.94825 -2 [-6, 2]; 0.14825
women 8 [6, 11]; >.999 11 [7, 14]; >.999 6 [2, 10]; 0.99875 -5 [-10, 1]; 0.04575
total 6 [4, 8]; >.999 8 [5, 10]; >.999 4 [2, 7]; >.999 -3 [-7, -0]; 0.0235
Sex difference 5 [2, 8]; 0.99825 6 [1, 11]; 0.994625 4 [-1, 9]; 0.93075 -2 [-9, 4]; 0.242125

8.7 Buffering of positive mood

f = "brmsP1/buffer/posaffect.Rdata"
the_data = my_data.Q_buffer[item_name %in% posaffect_items]
if (file.exists(f)) {
  load(f)
} else {
  fit = fit_ordinal_model(
    the_data, analysis_formula, f, family = acat())
}
Click to show sample diagnostics
Table 8.11: Summary of model paramters
parameter mean sd q5 q95 rhat ess_bulk ess_tail
Threshold[1] 0.2 0.59 -0.8 1.2 1.007 642 1348
Threshold[2] -1.5 0.58 -2.5 -0.6 1.007 622 1242
Threshold[3] -0.9 0.56 -1.8 0.0 1.007 570 1239
Threshold[4] -0.5 0.55 -1.4 0.4 1.007 572 1077
Threshold[5] -0.5 0.56 -1.4 0.4 1.008 567 1188
Threshold[6] -0.5 0.55 -1.4 0.4 1.008 584 1126
Threshold[7] -0.8 0.55 -1.7 0.1 1.008 541 1078
Threshold[8] -0.2 0.54 -1.1 0.7 1.008 548 1089
Threshold[9] -1.2 0.53 -2.0 -0.3 1.008 536 1015
Threshold[10] 0.2 0.53 -0.7 1.1 1.009 514 1047
Threshold[11] -0.2 0.53 -1.1 0.7 1.008 535 1038
Threshold[12] -0.2 0.53 -1.1 0.7 1.009 519 1020
Threshold[13] -0.1 0.53 -1.0 0.7 1.009 520 975
Threshold[14] 0.1 0.53 -0.7 1.0 1.009 519 1037
Threshold[15] -0.1 0.53 -0.9 0.8 1.008 518 963
Threshold[16] 0.4 0.53 -0.5 1.3 1.008 512 1041
Threshold[17] 0.2 0.53 -0.6 1.1 1.009 517 1031
Threshold[18] 0.6 0.53 -0.2 1.5 1.008 530 1020
Threshold[19] -0.5 0.53 -1.3 0.4 1.009 516 1028
b_Sexman:Drugoxycodone:Stresscontrol:reminderpre 0.2 0.52 -0.6 1.1 1.009 473 967
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpre -0.1 0.52 -0.9 0.8 1.007 472 1013
b_Sexman:Drugplacebo:Stresscontrol:reminderpre 0.1 0.51 -0.7 1.0 1.009 475 978
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpre 0.0 0.51 -0.8 0.9 1.007 475 1034
b_Sexman:Drugoxycodone:Stressstress:reminderpre 0.1 0.52 -0.7 1.0 1.009 477 966
b_Sexwoman:Drugoxycodone:Stressstress:reminderpre -0.2 0.51 -1.1 0.6 1.008 473 1035
b_Sexman:Drugplacebo:Stressstress:reminderpre 0.1 0.51 -0.7 1.0 1.009 476 1002
b_Sexwoman:Drugplacebo:Stressstress:reminderpre -0.1 0.51 -1.0 0.8 1.007 474 1034
b_Sexman:Drugoxycodone:Stresscontrol:reminderpost 0.2 0.52 -0.6 1.1 1.009 471 972
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpost 0.0 0.52 -0.8 0.9 1.007 477 1029
b_Sexman:Drugplacebo:Stresscontrol:reminderpost 0.1 0.51 -0.7 1.0 1.009 473 989
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpost 0.0 0.51 -0.8 0.9 1.008 475 1041
b_Sexman:Drugoxycodone:Stressstress:reminderpost 0.0 0.51 -0.8 0.9 1.009 477 1019
b_Sexwoman:Drugoxycodone:Stressstress:reminderpost -0.3 0.52 -1.1 0.6 1.008 475 1035
b_Sexman:Drugplacebo:Stressstress:reminderpost 0.0 0.51 -0.8 0.9 1.009 475 979
b_Sexwoman:Drugplacebo:Stressstress:reminderpost -0.3 0.51 -1.1 0.6 1.007 475 1044
sd_item_name__Intercept 0.3 0.14 0.1 0.5 1.001 1752 2514
sd_participant__Intercept 0.3 0.04 0.2 0.4 1.003 901 3238
sd_participant__Drugoxycodone:Stresscontrol:reminderpre 0.3 0.05 0.3 0.4 1.005 1005 1322
sd_participant__Drugplacebo:Stresscontrol:reminderpre 0.1 0.06 0.0 0.2 1.028 170 884
sd_participant__Drugoxycodone:Stressstress:reminderpre 0.3 0.05 0.2 0.4 1.010 463 850
sd_participant__Drugplacebo:Stressstress:reminderpre 0.1 0.05 0.0 0.2 1.016 262 1300
sd_participant__Drugoxycodone:Stresscontrol:reminderpost 0.3 0.05 0.2 0.4 1.004 1133 2156
sd_participant__Drugplacebo:Stresscontrol:reminderpost 0.1 0.06 0.0 0.2 1.028 181 1039
sd_participant__Drugoxycodone:Stressstress:reminderpost 0.3 0.05 0.2 0.4 1.012 362 836
sd_participant__Drugplacebo:Stressstress:reminderpost 0.2 0.05 0.1 0.3 1.022 216 1212
sd_session__Intercept 0.0 0.05 0.0 0.1 1.004 2208 3628
model: response.c ~ Sex:Drug:Stress:reminder + (1 | session) + (Drug:Stress:reminder | participant) + (1 | item_name)
b_ = fix effect parameters
sd_ standard deviation of hierarchical parameters
ess = effective sample size


p = plot_buffer_effect(fit,the_data,"posaffect")
((p$stress.plot + 
    ylab("stress effect on negative mood ratings")) | p$buffer.plot) + 
  plot_annotation(tag_levels = "A") + 
  plot_layout(widths = c(2,3))
A: Model-predicted effects of the first stress reminder on 'negative mood', 'selfconscious' and 'embarrased' ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Posterior distribution of the interaction of stress effect (pre-post reminder) ratings with drug condition. Dashed vertical lines enclode a region of practical equivalence.

Figure 8.12: A: Model-predicted effects of the first stress reminder on ‘negative mood’, ‘selfconscious’ and ‘embarrased’ ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Posterior distribution of the interaction of stress effect (pre-post reminder) ratings with drug condition. Dashed vertical lines enclode a region of practical equivalence.

p$buffer.stats %>% 
  kable(caption = "Buffering effect of oxycodone for negative mood ratings after reminders.") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic()
Table 8.12: Buffering effect of oxycodone for negative mood ratings after reminders.
Sex stress stress in placebo stress in oxycodone buffer
men -4 [-7, -1]; 0.0075 -6 [-10, -1]; 0.00775 -2 [-7, 2]; 0.163 3 [-3, 10]; 0.847375
women -8 [-12, -5]; <.001 -10 [-15, -5]; <.001 -7 [-12, -2]; 0.001875 3 [-5, 9]; 0.761375
total -6 [-9, -4]; <.001 -8 [-11, -4]; <.001 -5 [-8, -1]; 0.003 3 [-2, 8]; 0.883375
Sex difference -5 [-9, 0]; 0.0325 -4 [-11, 3]; 0.119625 -5 [-11, 2]; 0.0775 -1 [-10, 8]; 0.437125

9 Drug effects in Q5 and Q6

We start again with a look at the data

my_data.Q.drug56 =
  my_data.Q[,.(participant, session, stage, response, state.condition, Sex, Drug, item_name)] %>% 
  .[, session := factor(session)]

stage.names.b = gsub("\\\n"," ", stage.names)
my_data.Q.drug56 %>% 
  .[item_name %in% c("drug_like","drug_again","drug_dislike","high")] %>% 
  .[,.(m = mean(response), se = sd(response)/sqrt(length(unique(participant)))), 
    by = .(stage,Sex,state.condition, Drug, item_name)] %>% 
  .[, `:=`(lower = m-2*se, upper = m+2*se)] %>% 
  ggplot(aes(x = stage, y = m, color = Drug, lty = Sex, fill = Drug)) + 
  geom_line() + 
  facet_grid(item_name~state.condition, scale = "free") + 
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .2, color = NA) + 
  scale_col_drug + 
  scale_x_continuous(breaks = 1:11, labels = stage.names.b) + 
  theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust=1))

9.1 Drug liking

Here is a plot of drug liking, measured as the difference between the oxycodone and placebo conditions, dependent on sex.

my_data.Q.drug56 = 
my_data.Q.drug56[stage %in% 5:6]
the_data = my_data.Q.drug56[item_name == "drug_like"]
the_data %>% plot_drug_lines()
Drug liking stratified by stress condition and sex. Drug liking is measured as the slope of the lines from placebo to ocxcodone. Each line represents one participant.

Figure 9.1: Drug liking stratified by stress condition and sex. Drug liking is measured as the slope of the lines from placebo to ocxcodone. Each line represents one participant.

We can again summarize the data by plotting the distribution of the difference in liking between drug conditions.

p = plot_drug(the_data)
the_data = the_data %>% coarsen(var = "response", bins = 20, range = c(0,100))
p
Distribution of drug liking effects, stratified by stress condition (in panles) and sex (in colors).

Figure 9.2: Distribution of drug liking effects, stratified by stress condition (in panles) and sex (in colors).

raw_contrasts.drug = p$data[, .(participant, Sex, drug_effect, state.condition)][, Outcome := "drug_like"]

Here is the full factorial model we use to estimate the effect of stress:

drug_formula = response.c ~ state.condition:Sex:Drug + (1|session) + (1|participant/stage) + (0+state.condition:Drug|participant)

We estimate the model

fn = "brmsP1/drug/like.Rdata"
if (file.exists(fn)) {
  load(fn)
} else {
  fit = fit_ordinal_model(the_data, drug_formula, fn, family = acat())
}
#pp_check(fit, ndraws = 2000, type = "bars")
Click to show sample diagnostics
Table 9.1: Summary of model paramters
parameter mean sd q5 q95 rhat ess_bulk ess_tail
Threshold[1] 0.2 0.97 -1.5 1.8 1.002 772 1375
Threshold[2] 0.3 0.97 -1.3 1.9 1.001 894 1344
Threshold[3] -0.5 0.99 -2.1 1.1 1.001 878 1497
Threshold[4] -0.2 0.98 -1.8 1.5 1.002 882 1656
Threshold[5] 0.2 0.98 -1.5 1.8 1.001 887 1451
Threshold[6] 0.3 0.99 -1.3 1.9 1.001 946 1600
Threshold[7] 1.1 1.05 -0.6 2.8 1.001 1089 1765
Threshold[8] 0.0 1.07 -1.8 1.8 1.002 1058 1824
Threshold[9] 0.3 1.04 -1.4 2.0 1.001 994 1655
Threshold[10] 1.1 1.06 -0.7 2.9 1.001 1042 1970
Threshold[11] -0.2 1.03 -1.9 1.5 1.001 1056 1724
Threshold[12] 1.6 1.05 -0.1 3.3 1.001 1046 1849
Threshold[13] 0.2 1.04 -1.5 1.9 1.001 1089 1884
Threshold[14] 0.9 1.00 -0.8 2.6 1.002 960 1625
Threshold[15] 1.7 1.00 0.1 3.3 1.001 1034 1733
Threshold[16] 1.6 1.03 -0.1 3.3 1.002 967 1895
Threshold[17] 2.6 1.12 0.8 4.5 1.002 1267 2048
Threshold[18] 1.2 1.14 -0.7 3.1 1.001 1317 2535
Threshold[19] 1.0 1.02 -0.7 2.7 1.002 1042 1669
b_state.conditioncontrol:Sexman:Drugplacebo -1.6 0.89 -3.1 -0.2 1.002 1400 2190
b_state.conditionstress:Sexman:Drugplacebo -1.0 0.86 -2.5 0.3 1.000 1448 2054
b_state.conditioncontrol:Sexwoman:Drugplacebo -2.0 0.92 -3.6 -0.6 1.005 1350 2414
b_state.conditionstress:Sexwoman:Drugplacebo -1.8 0.94 -3.4 -0.3 1.002 1695 2252
b_state.conditioncontrol:Sexman:Drugoxycodone 2.0 0.78 0.8 3.3 1.000 1165 1995
b_state.conditionstress:Sexman:Drugoxycodone 1.8 0.79 0.5 3.1 1.000 1249 1781
b_state.conditioncontrol:Sexwoman:Drugoxycodone 1.3 0.76 0.0 2.6 1.001 1153 1611
b_state.conditionstress:Sexwoman:Drugoxycodone 1.4 0.79 0.1 2.7 1.001 1172 1970
sd_participant__Intercept 0.9 0.28 0.4 1.3 1.010 366 433
sd_participant__state.conditioncontrol:Drugplacebo 1.4 0.53 0.6 2.4 1.003 943 1029
sd_participant__state.conditionstress:Drugplacebo 1.7 0.49 1.0 2.6 1.006 982 1644
sd_participant__state.conditioncontrol:Drugoxycodone 0.8 0.35 0.2 1.4 1.009 323 542
sd_participant__state.conditionstress:Drugoxycodone 1.4 0.34 0.9 2.0 1.005 533 1029
sd_participant:stage__Intercept 0.1 0.07 0.0 0.2 1.002 1558 2235
sd_session__Intercept 1.8 0.87 0.5 3.3 1.003 642 712
model: response.c ~ state.condition:Sex:Drug + (1 | session) + (1 | participant/stage) + (0 + state.condition:Drug | participant)
b_ = fix effect parameters
sd_ standard deviation of hierarchical parameters
ess = effective sample size


For the following plots we calculate the effect of stress on drug liking as liking_in_oxycodone - liking_in_placebo, so that positive numbers mean that liking is greater in the stress condition.

Here are the estimated responses in the “cells” of the design, stratified by sex:

post.contrasts = get_drug_contrasts(fit, the_data)
attr(post.contrasts,"cell_stats") %>% 
  kable(caption = "Estimated responses for 'drug liking' by sex, state condition, drug") %>% 
  kable_styling(full_width = TRUE) %>% 
  kable_classic()
Table 9.2: Estimated responses for ‘drug liking’ by sex, state condition, drug
Sex state.condition placebo oxycodone
men control 2 (1, 3) 44 (41, 48)
men stress 7 (5, 9) 41 (38, 43)
women control 2 (1, 3) 27 (24, 30)
women stress 3 (2, 4) 34 (31, 36)
post.contrasts = get_drug_contrasts(fit, the_data)
drug_plotter(post.contrasts)
Drug liking. Left: Drug liking by sex. Top right: Drug liking by stress condition. Note that drug liking is the difference in liking between the ocycodone and placebo conditions.  Bottom right, clockwise from top left: Sex differences in drug liking, stress effect on drug liking, sex difference in drug liking in the stress condition, sex difference in drug liking in the control condition (for stress). Bottom right: Difference between the sex differences in liking between the control and stress condition.

Figure 9.3: Drug liking. Left: Drug liking by sex. Top right: Drug liking by stress condition. Note that drug liking is the difference in liking between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in drug liking, stress effect on drug liking, sex difference in drug liking in the stress condition, sex difference in drug liking in the control condition (for stress). Bottom right: Difference between the sex differences in liking between the control and stress condition.

drug_effects.stats = get_drug_effects.stats(post.contrasts, "drug_like")
drug_effects.stats.b = get_drug_effects.stats.b(post.contrasts, "drug_like")
make_tbl(post.contrasts, 
         caption = "Main and interaction effects of Drug on 'drug liking'.") %>%  
  add_footnote(c("Single words identify main effects in different conditions or sexes.\n '-' identifies contasts. 'in' signifies that an effect is reported for a subgroub or condition.","Numbers show drug effects."), notation="alphabet")
Table 9.3: Main and interaction effects of Drug on ‘drug liking’.
Effecta mean (CrI)b
women in control 25 [22, 28]; >.999
men in control 42 [39, 46]; >.999
women in stress 31 [28, 34]; >.999
men in stress 34 [31, 37]; >.999
women 28 [26, 30]; >.999
men 38 [36, 40]; >.999
stress 32 [30, 34]; >.999
control 34 [31, 36]; >.999
women - men -10 [-14, -7]; <.001
stress - ctrl -1 [-4, 2]; 0.2225
stress - ctrl in men -9 [-13, -4]; <.001
stress - ctrl in women 6 [2, 11]; 0.99725
women - men in ctrl -18 [-23, -12]; <.001
women - men in stress -3 [-7, 2]; 0.1165
three_way -15 [-22, -8]; <.001
a Single words identify main effects in different conditions or sexes.
‘-’ identifies contasts. ‘in’ signifies that an effect is reported for a subgroub or condition.
b Numbers show drug effects.

9.2 Drug Again

Raw data for the effect of drug.

Figure 9.4: Raw data for the effect of drug.

Table 9.4: Estimated responses for drug_again by sex, state condition, drug
Sex state.condition placebo oxycodone
men control 3 (2, 5) 31 (26, 35)
men stress 10 (7, 13) 31 (27, 36)
women control 6 (4, 9) 27 (22, 33)
women stress 9 (6, 13) 28 (23, 34)
Posterior distributions for effect of drug.

Figure 9.5: Posterior distributions for effect of drug.

The figure shows Left: drug again by sex. Top right: drug again by stress condition. Note that drug again is the difference in drug again between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in drug again, stress effect on drug again, sex difference in drug again in the stress condition, sex difference in drug again in the control condition (for stress). Bottom right: Difference between the sex differences in drug again between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘drug_again’.
Effecta mean (CrI)b
women in control 21 [15, 27]; >.999
men in control 27 [22, 32]; >.999
women in stress 19 [12, 26]; >.999
men in stress 22 [16, 27]; >.999
women 20 [16, 25]; >.999
men 24 [21, 28]; >.999
stress 20 [16, 25]; >.999
control 24 [20, 28]; >.999
women - men -4 [-10, 1]; 0.062
stress - ctrl -4 [-9, 2]; 0.0995
stress - ctrl in men -6 [-13, 1]; 0.053
stress - ctrl in women -2 [-11, 8]; 0.368
women - men in ctrl -6 [-14, 1]; 0.0505
women - men in stress -2 [-11, 6]; 0.304
three_way -4 [-16, 7]; 0.2535

9.3 Drug Dislike

Raw data for the effect of drug.

Figure 9.6: Raw data for the effect of drug.

Table 9.4: Estimated responses for drug_dislike by sex, state condition, drug
Sex state.condition placebo oxycodone
men control 5 (3, 8) 11 (7, 15)
men stress 8 (5, 12) 16 (11, 21)
women control 3 (2, 5) 24 (19, 29)
women stress 5 (3, 8) 32 (26, 37)
Posterior distributions for effect of drug.

Figure 9.7: Posterior distributions for effect of drug.

The figure shows Left: drug dislike by sex. Top right: drug dislike by stress condition. Note that drug dislike is the difference in drug dislike between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in drug dislike, stress effect on drug dislike, sex difference in drug dislike in the stress condition, sex difference in drug dislike in the control condition (for stress). Bottom right: Difference between the sex differences in drug dislike between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘drug_dislike’.
Effecta mean (CrI)b
women in control 21 [16, 26]; >.999
men in control 6 [1, 11]; 0.99475
women in stress 26 [20, 33]; >.999
men in stress 8 [2, 14]; 0.99525
women 24 [20, 28]; >.999
men 7 [3, 10]; >.999
stress 17 [13, 21]; >.999
control 13 [10, 17]; >.999
women - men 17 [11, 22]; >.999
stress - ctrl 4 [-2, 9]; 0.91975
stress - ctrl in men 2 [-6, 10]; 0.705
stress - ctrl in women 6 [-3, 14]; 0.9145
women - men in ctrl 15 [8, 22]; >.999
women - men in stress 19 [10, 27]; >.999
three_way -4 [-15, 8]; 0.25625

9.4 High

Raw data for the effect of drug.

Figure 9.8: Raw data for the effect of drug.

Table 9.4: Estimated responses for high by sex, state condition, drug
Sex state.condition placebo oxycodone
men control 3 (2, 5) 52 (47, 56)
men stress 10 (7, 12) 51 (46, 55)
women control 4 (3, 6) 59 (55, 64)
women stress 6 (4, 9) 61 (56, 65)
Posterior distributions for effect of drug.

Figure 9.9: Posterior distributions for effect of drug.

The figure shows Left: high by sex. Top right: high by stress condition. Note that high is the difference in high between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in high, stress effect on high, sex difference in high in the stress condition, sex difference in high in the control condition (for stress). Bottom right: Difference between the sex differences in high between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘high’.
Effecta mean (CrI)b
women in control 55 [50, 60]; >.999
men in control 49 [44, 53]; >.999
women in stress 54 [50, 59]; >.999
men in stress 41 [35, 46]; >.999
women 55 [51, 58]; >.999
men 45 [41, 48]; >.999
stress 48 [44, 51]; >.999
control 52 [48, 55]; >.999
women - men 10 [5, 15]; >.999
stress - ctrl -4 [-9, 1]; 0.0545
stress - ctrl in men -8 [-15, -1]; 0.01825
stress - ctrl in women -0 [-7, 6]; 0.44225
women - men in ctrl 6 [-1, 13]; 0.96575
women - men in stress 13 [6, 21]; >.999
three_way -7 [-17, 2]; 0.067

9.5 Dizzy

Raw data for the effect of drug.

Figure 9.10: Raw data for the effect of drug.

Table 9.4: Estimated responses for dizzy by sex, state condition, drug
Sex state.condition placebo oxycodone
men control 5 (3, 7) 40 (35, 46)
men stress 9 (6, 12) 33 (28, 38)
women control 8 (5, 11) 63 (58, 68)
women stress 11 (8, 14) 60 (54, 65)
Posterior distributions for effect of drug.

Figure 9.11: Posterior distributions for effect of drug.

The figure shows Left: dizzy by sex. Top right: dizzy by stress condition. Note that dizzy is the difference in dizzy between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in dizzy, stress effect on dizzy, sex difference in dizzy in the stress condition, sex difference in dizzy in the control condition (for stress). Bottom right: Difference between the sex differences in dizzy between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘dizzy’.
Effecta mean (CrI)b
women in control 55 [49, 61]; >.999
men in control 35 [30, 41]; >.999
women in stress 49 [43, 55]; >.999
men in stress 24 [18, 30]; >.999
women 52 [48, 56]; >.999
men 30 [26, 34]; >.999
stress 37 [32, 41]; >.999
control 45 [41, 49]; >.999
women - men 22 [16, 28]; >.999
stress - ctrl -9 [-15, -3]; 0.0025
stress - ctrl in men -11 [-20, -3]; 0.006
stress - ctrl in women -6 [-14, 3]; 0.082
women - men in ctrl 20 [11, 28]; >.999
women - men in stress 25 [16, 33]; >.999
three_way -5 [-16, 7]; 0.20875

9.6 Dry Mouth

Raw data for the effect of drug.

Figure 9.12: Raw data for the effect of drug.

Table 9.4: Estimated responses for dry_mouth by sex, state condition, drug
Sex state.condition placebo oxycodone
men control 7 (5, 10) 9 (6, 11)
men stress 14 (11, 18) 16 (13, 19)
women control 15 (12, 18) 14 (11, 17)
women stress 20 (16, 23) 24 (20, 27)
Posterior distributions for effect of drug.

Figure 9.13: Posterior distributions for effect of drug.

The figure shows Left: dry mouth by sex. Top right: dry mouth by stress condition. Note that dry mouth is the difference in dry mouth between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in dry mouth, stress effect on dry mouth, sex difference in dry mouth in the stress condition, sex difference in dry mouth in the control condition (for stress). Bottom right: Difference between the sex differences in dry mouth between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘dry_mouth’.
Effecta mean (CrI)b
women in control -1 [-5, 4]; 0.35
men in control 1 [-2, 4]; 0.75675
women in stress 4 [-1, 9]; 0.95025
men in stress 1 [-3, 6]; 0.73775
women 2 [-2, 5]; 0.849
men 1 [-2, 4]; 0.82425
stress 3 [-1, 6]; 0.947
control 0 [-3, 3]; 0.55075
women - men 0 [-4, 5]; 0.5595
stress - ctrl 3 [-2, 7]; 0.88225
stress - ctrl in men 0 [-5, 6]; 0.54675
stress - ctrl in women 5 [-2, 11]; 0.93125
women - men in ctrl -2 [-7, 3]; 0.2425
women - men in stress 3 [-4, 9]; 0.783
three_way -5 [-13, 4]; 0.14425

9.7 Blunted

Raw data for the effect of drug.

Figure 9.14: Raw data for the effect of drug.

Table 9.4: Estimated responses for blunted by sex, state condition, drug
Sex state.condition placebo oxycodone
men control 9 [7, 12]; 45 [40, 50];
men stress 14 [10, 17]; 35 [30, 40];
women control 15 [11, 18]; 57 [52, 62];
women stress 14 [11, 17]; 53 [48, 58];
Posterior distributions for effect of drug.

Figure 9.15: Posterior distributions for effect of drug.

The figure shows Left: blunted by sex. Top right: blunted by stress condition. Note that blunted is the difference in blunted between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in blunted, stress effect on blunted, sex difference in blunted in the stress condition, sex difference in blunted in the control condition (for stress). Bottom right: Difference between the sex differences in blunted between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘blunted’.
Effecta mean (CrI)b
women in control 42 [36, 49]; >.999
men in control 36 [30, 42]; >.999
women in stress 39 [33, 45]; >.999
men in stress 22 [16, 27]; >.999
women 41 [37, 45]; >.999
men 29 [25, 33]; >.999
stress 30 [26, 35]; >.999
control 39 [35, 43]; >.999
women - men 12 [6, 18]; >.999
stress - ctrl -9 [-15, -3]; 0.00325
stress - ctrl in men -15 [-23, -6]; <.001
stress - ctrl in women -3 [-12, 6]; 0.2495
women - men in ctrl 6 [-2, 15]; 0.92675
women - men in stress 18 [9, 26]; >.999
three_way -12 [-24, 0]; 0.02675

9.8 Nauseous

Raw data for the effect of drug.

Figure 9.16: Raw data for the effect of drug.

Table 9.4: Estimated responses for nauseous by sex, state condition, drug
Sex state.condition placebo oxycodone
men control 3 (2, 4) 9 (6, 13)
men stress 5 (3, 7) 6 (4, 9)
women control 3 (2, 4) 15 (11, 19)
women stress 4 (3, 5) 13 (10, 17)
Posterior distributions for effect of drug.

Figure 9.17: Posterior distributions for effect of drug.

The figure shows Left: nauseous by sex. Top right: nauseous by stress condition. Note that nauseous is the difference in nauseous between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in nauseous, stress effect on nauseous, sex difference in nauseous in the stress condition, sex difference in nauseous in the control condition (for stress). Bottom right: Difference between the sex differences in nauseous between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘nauseous’.
Effecta mean (CrI)b
women in control 12 [8, 16]; >.999
men in control 6 [3, 10]; >.999
women in stress 9 [6, 13]; >.999
men in stress 2 [-2, 5]; 0.823
women 11 [8, 14]; >.999
men 4 [2, 6]; >.999
stress 5 [3, 8]; >.999
control 9 [7, 12]; >.999
women - men 7 [3, 11]; >.999
stress - ctrl -4 [-7, -0]; 0.0225
stress - ctrl in men -5 [-10, -0]; 0.0215
stress - ctrl in women -3 [-9, 3]; 0.1635
women - men in ctrl 6 [1, 11]; 0.98575
women - men in stress 8 [3, 13]; 0.998
three_way -2 [-10, 6]; 0.32225

9.9 Warm Face

Raw data for the effect of drug.

Figure 9.18: Raw data for the effect of drug.

Table 9.4: Estimated responses for warm_face by sex, state condition, drug
Sex state.condition placebo oxycodone
men control 5 (3, 7) 24 (19, 28)
men stress 6 (4, 9) 21 (17, 25)
women control 5 (3, 6) 25 (20, 29)
women stress 11 (8, 14) 35 (30, 40)
Posterior distributions for effect of drug.

Figure 9.19: Posterior distributions for effect of drug.

The figure shows Left: warm face by sex. Top right: warm face by stress condition. Note that warm face is the difference in warm face between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in warm face, stress effect on warm face, sex difference in warm face in the stress condition, sex difference in warm face in the control condition (for stress). Bottom right: Difference between the sex differences in warm face between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘warm_face’.
Effecta mean (CrI)b
women in control 20 [16, 25]; >.999
men in control 19 [14, 24]; >.999
women in stress 24 [18, 29]; >.999
men in stress 15 [10, 19]; >.999
women 22 [18, 26]; >.999
men 17 [13, 20]; >.999
stress 19 [16, 23]; >.999
control 20 [16, 23]; >.999
women - men 5 [0, 10]; 0.98075
stress - ctrl -0 [-6, 5]; 0.43375
stress - ctrl in men -5 [-11, 2]; 0.0845
stress - ctrl in women 4 [-4, 11]; 0.8295
women - men in ctrl 1 [-6, 8]; 0.6125
women - men in stress 9 [2, 17]; 0.994
three_way -8 [-18, 2]; 0.05175

9.10 Not Myself

Raw data for the effect of drug.

Figure 9.20: Raw data for the effect of drug.

Table 9.4: Estimated responses for not_myself by sex, state condition, drug
Sex state.condition placebo oxycodone
men control 8 (5, 10) 28 (23, 33)
men stress 10 (8, 14) 31 (27, 36)
women control 8 (6, 11) 45 (40, 50)
women stress 10 (7, 14) 43 (38, 49)
Posterior distributions for effect of drug.

Figure 9.21: Posterior distributions for effect of drug.

The figure shows Left: not myself by sex. Top right: not myself by stress condition. Note that not myself is the difference in not myself between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in not myself, stress effect on not myself, sex difference in not myself in the stress condition, sex difference in not myself in the control condition (for stress). Bottom right: Difference between the sex differences in not myself between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘not_myself’.
Effecta mean (CrI)b
women in control 37 [31, 42]; >.999
men in control 21 [15, 26]; >.999
women in stress 33 [27, 39]; >.999
men in stress 21 [15, 27]; >.999
women 35 [31, 39]; >.999
men 21 [17, 25]; >.999
stress 27 [23, 31]; >.999
control 29 [25, 33]; >.999
women - men 14 [8, 20]; >.999
stress - ctrl -2 [-7, 4]; 0.28475
stress - ctrl in men 0 [-7, 8]; 0.53875
stress - ctrl in women -4 [-12, 5]; 0.19475
women - men in ctrl 16 [8, 24]; >.999
women - men in stress 12 [4, 20]; 0.99825
three_way 4 [-7, 15]; 0.75825

9.11 Pos Affect

Raw data for the effect of drug.

Figure 9.22: Raw data for the effect of drug.

Table 9.4: Estimated responses for pos_affect by sex, state condition, drug
Sex state.condition placebo oxycodone
men control 72 (70, 74) 73 (71, 75)
men stress 68 (66, 70) 69 (67, 71)
women control 68 (66, 70) 64 (62, 66)
women stress 57 (55, 59) 53 (51, 55)
Posterior distributions for effect of drug.

Figure 9.23: Posterior distributions for effect of drug.

The figure shows Left: pos affect by sex. Top right: pos affect by stress condition. Note that pos affect is the difference in pos affect between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in pos affect, stress effect on pos affect, sex difference in pos affect in the stress condition, sex difference in pos affect in the control condition (for stress). Bottom right: Difference between the sex differences in pos affect between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘pos_affect’.
Effecta mean (CrI)b
women in control -4 [-7, -1]; 0.00175
men in control 1 [-2, 3]; 0.66475
women in stress -4 [-7, -1]; 0.0055
men in stress 1 [-2, 4]; 0.789
women -4 [-6, -2]; <.001
men 1 [-1, 3]; 0.81025
stress -1 [-3, 1]; 0.08075
control -2 [-4, 0]; 0.03775
women - men -5 [-8, -2]; <.001
stress - ctrl 0 [-2, 3]; 0.57825
stress - ctrl in men 1 [-3, 4]; 0.60825
stress - ctrl in women 0 [-4, 4]; 0.508
women - men in ctrl -5 [-8, -1]; 0.0085
women - men in stress -5 [-9, -1]; 0.00575
three_way 0 [-5, 6]; 0.56825

9.12 Neg Affect

Raw data for the effect of drug.

Figure 9.24: Raw data for the effect of drug.

Table 9.4: Estimated responses for neg_affect by sex, state condition, drug
Sex state.condition placebo oxycodone
men control 4 (3, 5) 5 (4, 6)
men stress 10 (8, 11) 9 (7, 10)
women control 5 (4, 6) 9 (8, 10)
women stress 19 (17, 22) 18 (16, 20)
Posterior distributions for effect of drug.

Figure 9.25: Posterior distributions for effect of drug.

The figure shows Left: neg affect by sex. Top right: neg affect by stress condition. Note that neg affect is the difference in neg affect between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in neg affect, stress effect on neg affect, sex difference in neg affect in the stress condition, sex difference in neg affect in the control condition (for stress). Bottom right: Difference between the sex differences in neg affect between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘neg_affect’.
Effecta mean (CrI)b
women in control 4 [2, 6]; >.999
men in control 1 [-0, 2]; 0.962
women in stress -1 [-4, 2]; 0.1995
men in stress -1 [-3, 1]; 0.1375
women 1 [-0, 3]; 0.94125
men 0 [-1, 1]; 0.51775
stress -1 [-3, 1]; 0.1005
control 3 [2, 4]; >.999
women - men 1 [-1, 3]; 0.90625
stress - ctrl -4 [-6, -2]; <.001
stress - ctrl in men -2 [-5, 0]; 0.031
stress - ctrl in women -5 [-8, -2]; <.001
women - men in ctrl 3 [1, 5]; 0.99625
women - men in stress -0 [-4, 3]; 0.45525
three_way 3 [-1, 7]; 0.93375

9.13 Summary of Drug effects

drug_effects.stats %>% 
  .[Outcome == "high", Outcome := "feeling_high"] %>% 
  .[grepl("high|drug|dizzy|blunted|not|warm", Outcome)] %>% 
  .[,Outcome := stringr::str_to_title(Outcome)] %>% 
  .[, Outcome := gsub("_"," ", Outcome)] %>% 
  kable(caption = "Drug effects on subjective ratings.") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic() 
Table 9.5: Drug effects on subjective ratings.
Outcome Men Women Total Sex difference
Drug like 38 [36, 40]; 28 [26, 30]; 33 [31, 34]; -10 [-14, -7]; <.001
Drug again 24 [21, 28]; 20 [16, 25]; 22 [19, 25]; -4 [-10, 1]; 0.062
Drug dislike 7 [3, 10]; 24 [20, 28]; 15 [13, 18]; 17 [11, 22]; >.999
Feeling high 45 [41, 48]; 55 [51, 58]; 50 [47, 52]; 10 [5, 15]; >.999
Dizzy 30 [26, 34]; 52 [48, 56]; 41 [38, 44]; 22 [16, 28]; >.999
Blunted 29 [25, 33]; 41 [37, 45]; 35 [32, 38]; 12 [6, 18]; >.999
Warm face 17 [13, 20]; 22 [18, 26]; 19 [17, 22]; 5 [0, 10]; 0.98075
Not myself 21 [17, 25]; 35 [31, 39]; 28 [25, 31]; 14 [8, 20]; >.999
clean = function(dt) {
  dt = 
    dt %>%  
    .[grepl("lik|high", Outcome)] %>%  
    .[Outcome == "high", Outcome := "feeling_high"] %>% 
    .[,Outcome := stringr::str_to_title(Outcome)] %>%
    .[,Sex := ifelse(Sex == "men","Men","Women")] %>%
    .[,state.condition := stringr::str_to_title(state.condition)] %>%
    .[, Outcome := gsub("_"," ", Outcome)] %>% 
    .[, Outcome := gsub("like","liking", Outcome)] %>% 
    .[, Outcome := paste("\u0394 ",Outcome)] %>% 
    .[, Outcome := factor(Outcome)]
  o.levels = levels(dt$Outcome)[c(3,2,1)]
  dt[, Outcome := factor(Outcome, levels = o.levels)]
  return(dt)
}

raw_contrasts.drug.p = 
  raw_contrasts.drug %>%
  copy() %>% 
  clean() %>%  
  .[, l := fmean(drug_effect)-3*fsd(drug_effect), by = .(Outcome,Sex)] %>% 
  .[, u := fmean(drug_effect)+3*fsd(drug_effect), by = .(Outcome,Sex)] %>% 
  .[, plot := drug_effect < u & drug_effect > l]

breaks_drug_fun <- function(x) {
  d.count <<- d.count + 1L
  switch(
    d.count,
    seq(0,100,25), c(0),
    seq(-25,100,25),c(0),
    seq(-50,100,50), c(0),
    seq(0,100,25), c(0),
    seq(-25,100,25),c(0),
    seq(-50,100,50),
  )
}

d.count = 0
# summary_plot.drug = 
#   drug_effects.stats.b %>% 
#     clean() %>% 
#     ggplot(aes(x = state.condition, y = drug_effect,
#                fill = Sex, shape = state.condition,
#                group = factor(Sex):factor(state.condition))) + 
#     geom_bar(stat = "identity", position = position_dodge(1)) + 
#     geom_hline(yintercept = 0) +
#     geom_errorbar(aes(ymin = lower, ymax = upper),  linewidth = 1,
#                   width = 1/5, position = position_dodge(1)) + 
#     geom_quasirandom(data = raw_contrasts.drug.p[plot == TRUE],
#                      dodge.width = 1,
#                      color = adjustcolor("black",alpha = .5),
#                      aes(shape = state.condition)) +
#     scale_shape_manual(values = c(1,2)) + 
#     facet_grid(Outcome ~ Sex, scale = "free_y",switch = "y") + 
#     #scale_y_continuous(breaks = breaks_drug_fun()) + 
#     guides(fill = "none", 
#            shape = guide_legend(nrow = 1, title.position = "left", title = "State induction")) + 
#     theme(legend.position = "top", axis.line.y = element_line(),
#           panel.grid.major.y = element_blank(), panel.spacing = unit(1, "lines"),
#           strip.background = element_blank(), legend.key.size = unit(1,"line"),
#           strip.placement = "outside", legend.title=element_text(size=12)) + 
#     scale_col_sex2 + ylab("") + xlab("")
summary_plot.drug = 
  drug_effects.stats.b %>% 
  clean() %>% 
  ggplot(aes(x = Sex, y = drug_effect,
             fill = Sex, shape = state.condition,
             group = factor(Sex):factor(state.condition))) + 
  geom_bar(stat = "identity", position = position_dodge(1)) + 
  geom_hline(yintercept = 0) +
  geom_errorbar(aes(ymin = lower, ymax = upper),  linewidth = 1,
                width = 1/5, position = position_dodge(1)) + 
  geom_quasirandom(data = raw_contrasts.drug.p[plot == TRUE],
                   dodge.width = 1,
                   color = adjustcolor("black",alpha = .5),
                   aes(shape = state.condition)) +
  scale_shape_manual(values = c(1,2)) + 
  facet_grid(Outcome ~ ., scale = "free_y",switch = "y") + 
  #scale_y_continuous(breaks = breaks_drug_fun()) + 
  guides(fill = "none", 
         shape = guide_legend(nrow = 1, title.position = "left", title = "State induction")) + 
  theme(legend.position = "top", axis.line.y = element_line(),
        panel.grid.major.y = element_blank(), panel.spacing = unit(1, "lines"),
        strip.background = element_blank(), legend.key.size = unit(1,"line"),
        strip.placement = "outside", legend.title=element_text(size=12)) + 
  scale_col_sex2 + ylab("") + xlab("") + gg_no_x_axis
emb_plot = 
  select_data("embarassed") %>% 
    plot_stages(
        my_stages = 2:9,
        by.Sex = TRUE, lty = "Sex",
        shape = "state.condition", 
        by.Drug = TRUE, fill = "Drug") + 
    facet_wrap(~Sex, nrow = 2) + 
    theme(strip.background = element_blank(),
          legend.position = c(.5,.9),
          legend.margin=margin(0,0,0,0),
          legend.spacing = unit(.25,"line"),
          legend.title=element_text(size=12)) + 
    coord_cartesian(ylim = c(0,100), xlim = c(-20,35)) + 
    guides(fill = guide_legend(nrow = 2, title.position = "top", title = "Drug"), 
           shape = "none",
           lty = "none", color = "none") +
    ylab("Embarrassed (VAS)")
summary_plot = 
  (summary_plot.stress | summary_plot.drug | emb_plot) + 
  plot_layout(widths = c(1,2,1.5)) +
  plot_annotation(tag_levels = "A")
s.count = 0
#d.count = 0
ggsave(summary_plot,file = "figures/summary_plot.png", width = 9, height = 7)
s.count = 0
#d.count = 0
ggsave(summary_plot,file = "figures/summary_plot.pdf", width = 9, height = 7)
Overview over effects of stress and drug
Overview over effects of stress and drug

  1. For instance, when examining the stress effect dependent of drug, we have 8 conditions (2 drug conditions time 2 stress conditions times 2 measurements pre and post stress induction). Then we calculate for each of these 8 conditions the response probability for each of the 20 response levels.↩︎